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

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

minor changes to implement Aboca

File size: 150.1 KB
Line 
1/* $Id: CoinAbcBaseFactorization4.cpp 1373 2011-01-03 23:57:44Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5#if CLP_HAS_ABC
6#ifdef ABC_JUST_ONE_FACTORIZATION
7#include "CoinAbcCommonFactorization.hpp"
8#define CoinAbcTypeFactorization CoinAbcBaseFactorization
9#define ABC_SMALL -1
10#include "CoinAbcBaseFactorization.hpp"
11#endif
12#ifdef CoinAbcTypeFactorization
13
14#include "CoinIndexedVector.hpp"
15#include "CoinHelperFunctions.hpp"
16#include "CoinAbcHelperFunctions.hpp"
17#if CILK_CONFLICT>0
18// for conflicts
19extern int cilk_conflict;
20#endif
21
22//:class CoinAbcTypeFactorization.  Deals with Factorization and Updates
23#if ABC_SMALL<2
24//  getColumnSpaceIterateR.  Gets space for one extra R element in Column
25//may have to do compression  (returns true)
26//also moves existing vector
27bool
28CoinAbcTypeFactorization::getColumnSpaceIterateR ( CoinSimplexInt iColumn, CoinFactorizationDouble value,
29                                           CoinSimplexInt iRow)
30{
31  CoinFactorizationDouble * COIN_RESTRICT elementR = elementRAddress_ + lengthAreaR_;
32  CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_ + lengthAreaR_;
33  CoinBigIndex * COIN_RESTRICT startR = startColumnRAddress_+maximumPivots_+1;
34  CoinSimplexInt * COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
35  CoinSimplexInt number = numberInColumnPlus[iColumn];
36  //*** modify so sees if can go in
37  //see if it can go in at end
38  CoinSimplexInt * COIN_RESTRICT nextColumn = nextColumnAddress_;
39  CoinSimplexInt * COIN_RESTRICT lastColumn = lastColumnAddress_;
40  if (lengthAreaR_-startR[maximumRowsExtra_]<number+1) {
41    //compression
42    CoinSimplexInt jColumn = nextColumn[maximumRowsExtra_];
43    CoinBigIndex put = 0;
44    while ( jColumn != maximumRowsExtra_ ) {
45      //move
46      CoinBigIndex get;
47      CoinBigIndex getEnd;
48      get = startR[jColumn];
49      getEnd = get + numberInColumnPlus[jColumn];
50      startR[jColumn] = put;
51      CoinBigIndex i;
52      for ( i = get; i < getEnd; i++ ) {
53        indexRowR[put] = indexRowR[i];
54        elementR[put] = elementR[i];
55        put++;
56      }
57      jColumn = nextColumn[jColumn];
58    }
59    numberCompressions_++;
60    startR[maximumRowsExtra_]=put;
61  }
62  // Still may not be room (as iColumn was still in)
63  if (lengthAreaR_-startR[maximumRowsExtra_]<number+1) 
64    return false;
65
66  CoinSimplexInt next = nextColumn[iColumn];
67  CoinSimplexInt last = lastColumn[iColumn];
68  //out
69  nextColumn[last] = next;
70  lastColumn[next] = last;
71
72  CoinBigIndex put = startR[maximumRowsExtra_];
73  //in at end
74  last = lastColumn[maximumRowsExtra_];
75  nextColumn[last] = iColumn;
76  lastColumn[maximumRowsExtra_] = iColumn;
77  lastColumn[iColumn] = last;
78  nextColumn[iColumn] = maximumRowsExtra_;
79 
80  //move
81  CoinBigIndex get = startR[iColumn];
82  startR[iColumn] = put;
83  CoinSimplexInt i = 0;
84  for (i=0 ; i < number; i ++ ) {
85    elementR[put]= elementR[get];
86    indexRowR[put++] = indexRowR[get++];
87  }
88  //insert
89  elementR[put]=value;
90  indexRowR[put++]=iRow;
91  numberInColumnPlus[iColumn]++;
92  //add 4 for luck
93  startR[maximumRowsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaR_);
94  return true;
95}
96#endif
97CoinSimplexInt CoinAbcTypeFactorization::checkPivot(CoinSimplexDouble saveFromU,
98                                 CoinSimplexDouble oldPivot) const
99{
100  CoinSimplexInt status;
101  if ( fabs ( saveFromU ) > 1.0e-8 ) {
102    CoinFactorizationDouble checkTolerance;
103    if ( numberRowsExtra_ < numberRows_ + 2 ) {
104      checkTolerance = 1.0e-5;
105    } else if ( numberRowsExtra_ < numberRows_ + 10 ) {
106      checkTolerance = 1.0e-6;
107    } else if ( numberRowsExtra_ < numberRows_ + 50 ) {
108      checkTolerance = 1.0e-8;
109    } else {
110      checkTolerance = 1.0e-10;
111    }       
112    checkTolerance *= relaxCheck_;
113    if ( fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < checkTolerance ) {
114      status = 0;
115    } else {
116#if COIN_DEBUG
117      std::cout <<"inaccurate pivot "<< oldPivot << " " 
118                << saveFromU << std::endl;
119#endif
120      if ( fabs ( fabs ( oldPivot ) - fabs ( saveFromU ) ) < 1.0e-12 ||
121        fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < 1.0e-8 ) {
122        status = 1;
123      } else {
124        status = 2;
125      }       
126    }       
127  } else {
128    //error
129    status = 2;
130#if COIN_DEBUG
131    std::cout <<"inaccurate pivot "<< saveFromU / oldPivot
132              << " " << saveFromU << std::endl;
133#endif
134  } 
135  return status;
136}
137#define INLINE_IT
138#define INLINE_IT2
139#define INLINE_IT3
140#define UNROLL 0
141inline void scatterUpdateInline(CoinSimplexInt number,
142                          CoinFactorizationDouble pivotValue,
143                          const CoinFactorizationDouble *  COIN_RESTRICT thisElement,
144                          const CoinSimplexInt *  COIN_RESTRICT thisIndex,
145                          CoinFactorizationDouble *  COIN_RESTRICT region)
146{
147#if UNROLL==0
148  for (CoinBigIndex j=number-1 ; j >=0; j-- ) {
149    CoinSimplexInt iRow = thisIndex[j];
150    CoinFactorizationDouble regionValue = region[iRow];
151    CoinFactorizationDouble value = thisElement[j];
152    assert (value);
153    region[iRow] = regionValue - value * pivotValue;
154  }
155#elif UNROLL==1
156  if ((number&1)!=0) {
157    number--;
158    CoinSimplexInt iRow = thisIndex[number];
159    CoinFactorizationDouble regionValue = region[iRow];
160    CoinFactorizationDouble value = thisElement[number];
161    region[iRow] = regionValue - value * pivotValue;
162  }
163  for (CoinBigIndex j=number-1 ; j >=0; j-=2 ) {
164    CoinSimplexInt iRow0 = thisIndex[j];
165    CoinSimplexInt iRow1 = thisIndex[j-1];
166    CoinFactorizationDouble regionValue0 = region[iRow0];
167    CoinFactorizationDouble regionValue1 = region[iRow1];
168    region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
169    region[iRow1] = regionValue1 - thisElement[j-1] * pivotValue;
170  }
171#elif UNROLL==2
172  CoinSimplexInt iRow0;
173  CoinSimplexInt iRow1;
174  CoinFactorizationDouble regionValue0;
175  CoinFactorizationDouble regionValue1;
176  switch(number) {
177  case 0:
178    break;
179  case 1:
180    iRow0 = thisIndex[0];
181    regionValue0 = region[iRow0];
182    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
183    break;
184  case 2:
185    iRow0 = thisIndex[0];
186    iRow1 = thisIndex[1];
187    regionValue0 = region[iRow0];
188    regionValue1 = region[iRow1];
189    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
190    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
191    break;
192  case 3:
193    iRow0 = thisIndex[0];
194    iRow1 = thisIndex[1];
195    regionValue0 = region[iRow0];
196    regionValue1 = region[iRow1];
197    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
198    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
199    iRow0 = thisIndex[2];
200    regionValue0 = region[iRow0];
201    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
202    break;
203  case 4:
204    iRow0 = thisIndex[0];
205    iRow1 = thisIndex[1];
206    regionValue0 = region[iRow0];
207    regionValue1 = region[iRow1];
208    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
209    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
210    iRow0 = thisIndex[2];
211    iRow1 = thisIndex[3];
212    regionValue0 = region[iRow0];
213    regionValue1 = region[iRow1];
214    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
215    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
216    break;
217  case 5:
218    iRow0 = thisIndex[0];
219    iRow1 = thisIndex[1];
220    regionValue0 = region[iRow0];
221    regionValue1 = region[iRow1];
222    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
223    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
224    iRow0 = thisIndex[2];
225    iRow1 = thisIndex[3];
226    regionValue0 = region[iRow0];
227    regionValue1 = region[iRow1];
228    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
229    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
230    iRow0 = thisIndex[4];
231    regionValue0 = region[iRow0];
232    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
233    break;
234  case 6:
235    iRow0 = thisIndex[0];
236    iRow1 = thisIndex[1];
237    regionValue0 = region[iRow0];
238    regionValue1 = region[iRow1];
239    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
240    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
241    iRow0 = thisIndex[2];
242    iRow1 = thisIndex[3];
243    regionValue0 = region[iRow0];
244    regionValue1 = region[iRow1];
245    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
246    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
247    iRow0 = thisIndex[4];
248    iRow1 = thisIndex[5];
249    regionValue0 = region[iRow0];
250    regionValue1 = region[iRow1];
251    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
252    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
253    break;
254  case 7:
255    iRow0 = thisIndex[0];
256    iRow1 = thisIndex[1];
257    regionValue0 = region[iRow0];
258    regionValue1 = region[iRow1];
259    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
260    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
261    iRow0 = thisIndex[2];
262    iRow1 = thisIndex[3];
263    regionValue0 = region[iRow0];
264    regionValue1 = region[iRow1];
265    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
266    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
267    iRow0 = thisIndex[4];
268    iRow1 = thisIndex[5];
269    regionValue0 = region[iRow0];
270    regionValue1 = region[iRow1];
271    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
272    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
273    iRow0 = thisIndex[6];
274    regionValue0 = region[iRow0];
275    region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
276    break;
277  case 8:
278    iRow0 = thisIndex[0];
279    iRow1 = thisIndex[1];
280    regionValue0 = region[iRow0];
281    regionValue1 = region[iRow1];
282    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
283    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
284    iRow0 = thisIndex[2];
285    iRow1 = thisIndex[3];
286    regionValue0 = region[iRow0];
287    regionValue1 = region[iRow1];
288    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
289    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
290    iRow0 = thisIndex[4];
291    iRow1 = thisIndex[5];
292    regionValue0 = region[iRow0];
293    regionValue1 = region[iRow1];
294    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
295    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
296    iRow0 = thisIndex[6];
297    iRow1 = thisIndex[7];
298    regionValue0 = region[iRow0];
299    regionValue1 = region[iRow1];
300    region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
301    region[iRow1] = regionValue1 - thisElement[7] * pivotValue;
302    break;
303  default:
304    if ((number&1)!=0) {
305      number--;
306      CoinSimplexInt iRow = thisIndex[number];
307      CoinFactorizationDouble regionValue = region[iRow];
308      CoinFactorizationDouble value = thisElement[number];
309      region[iRow] = regionValue - value * pivotValue;
310    }
311    for (CoinBigIndex j=number-1 ; j >=0; j-=2 ) {
312      CoinSimplexInt iRow0 = thisIndex[j];
313      CoinSimplexInt iRow1 = thisIndex[j-1];
314      CoinFactorizationDouble regionValue0 = region[iRow0];
315      CoinFactorizationDouble regionValue1 = region[iRow1];
316      region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
317      region[iRow1] = regionValue1 - thisElement[j-1] * pivotValue;
318    }
319    break;
320  }
321#endif
322}
323inline CoinFactorizationDouble gatherUpdate(CoinSimplexInt number,
324                          const CoinFactorizationDouble *  COIN_RESTRICT thisElement,
325                          const CoinSimplexInt *  COIN_RESTRICT thisIndex,
326                          CoinFactorizationDouble *  COIN_RESTRICT region)
327{
328  CoinFactorizationDouble pivotValue=0.0;
329  for (CoinBigIndex j = 0; j < number; j ++ ) {
330    CoinFactorizationDouble value = thisElement[j];
331    CoinSimplexInt jRow = thisIndex[j];
332    value *= region[jRow];
333    pivotValue -= value;
334  }
335  return pivotValue;
336}
337inline void multiplyIndexed(CoinSimplexInt number,
338                            const CoinFactorizationDouble *  COIN_RESTRICT multiplier,
339                            const CoinSimplexInt *  COIN_RESTRICT thisIndex,
340                            CoinFactorizationDouble *  COIN_RESTRICT region)
341{
342  for (CoinSimplexInt i = 0; i < number; i ++ ) {
343    CoinSimplexInt iRow = thisIndex[i];
344    CoinSimplexDouble value = region[iRow];
345    value *= multiplier[iRow];
346    region[iRow] = value;
347  }
348}
349#if 0
350/* Checks if can replace one Column to basis,
351   returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots
352   Fills in region for use later
353   partial update already in U */
354int
355CoinAbcTypeFactorization::checkReplace ( const AbcSimplex * model,
356                                         CoinIndexedVector * regionSparse,
357                                         int pivotRow,
358                                         CoinSimplexDouble &pivotCheck,
359                                         double acceptablePivot)
360{
361  if ( lengthR_+numberRows_ >= lengthAreaR_ ) {
362    //not enough room
363    return 3;
364  }       
365#ifdef ABC_USE_FUNCTION_POINTERS
366  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
367  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
368  CoinBigIndex startU = scatter[numberRows_].offset;
369  CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
370  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
371  CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
372#else
373  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
374  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
375  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
376  CoinBigIndex startU = startColumnU[numberRows_];
377  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
378  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
379  CoinSimplexInt numberInColumnU2 = numberInColumn[numberRows_];
380#endif
381  if ( lengthU_+numberInColumnU2 >= lengthAreaU_ ) {
382    //not enough room
383    return 3;
384  }
385#if ABC_SMALL<2
386  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
387#endif
388  //zeroed out region
389  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
390  CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
391  CoinFactorizationDouble oldPivot = pivotRegion[pivotRow];
392  // for accuracy check
393  pivotCheck = pivotCheck / oldPivot;
394
395#ifndef ABC_USE_FUNCTION_POINTERS
396  CoinBigIndex saveStart = startColumnU[pivotRow];
397  CoinBigIndex saveEnd = saveStart + numberInColumn[pivotRow];
398#endif
399  //get entries in row (pivot not stored)
400  CoinSimplexInt numberNonZero = 0;
401#if ABC_SMALL<2
402  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
403#endif
404#if CONVERTROW
405  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
406#if CONVERTROW>2
407  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
408#endif
409#endif
410#if ABC_SMALL<2
411  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
412  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
413#endif
414  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
415  CoinSimplexInt status=0;
416#ifndef ABC_USE_FUNCTION_POINTERS
417  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
418#endif
419
420#ifndef NDEBUG
421#if CONVERTROW
422#define CONVERTDEBUG
423#endif
424#endif
425  int nInRow;
426#if ABC_SMALL<2
427  CoinBigIndex start=0;
428  CoinBigIndex end=0;
429  if (gotUCopy()) {
430    start=startRowU[pivotRow];
431    nInRow=numberInRowU[pivotRow];
432    end= start + nInRow;
433    CoinSimplexInt smallestIndex=pivotLinkedForwardsAddress_[-1];
434    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
435    if (nInRow<10) {
436      CoinSimplexInt smallest=numberRowsExtra_;
437      for (CoinBigIndex i = start; i < end ; i ++ ) {
438        CoinSimplexInt jColumn = indexColumnU[i];
439        if (pivotRowForward[jColumn]<smallest) {
440          smallest=pivotRowForward[jColumn];
441          smallestIndex=jColumn;
442        }
443#ifndef ABC_USE_FUNCTION_POINTERS
444#ifdef CONVERTDEBUG
445        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
446        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
447#endif     
448        region[jColumn] = elementRowU[i];
449#else
450#ifdef CONVERTDEBUG
451        CoinBigIndex j = convertRowToColumn[i];
452        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
453        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
454#endif     
455        region[jColumn] = elementRowU[i];
456#endif
457        regionIndex[numberNonZero++] = jColumn;
458      }
459    } else {
460      for (CoinBigIndex i = start; i < end ; i ++ ) {
461        CoinSimplexInt jColumn = indexColumnU[i];
462#ifdef CONVERTDEBUG
463#ifdef ABC_USE_FUNCTION_POINTERS
464        CoinBigIndex j = convertRowToColumn[i];
465        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
466        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
467#else
468        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
469        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
470#endif
471#endif     
472        region[jColumn] = elementRowU[i];
473        regionIndex[numberNonZero++] = jColumn;
474      }
475    }
476    //do BTRAN - finding first one to use
477    regionSparse->setNumElements ( numberNonZero );
478    if (numberNonZero) {
479      assert (smallestIndex>=0);
480#ifndef NDEBUG
481      const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
482      CoinSimplexInt jRow=pivotLinked[-1];
483      if (jRow!=smallestIndex) {
484        while (jRow>=0) {
485          CoinFactorizationDouble pivotValue = region[jRow];
486          if (pivotValue)
487            break;
488          jRow=pivotLinked[jRow];
489        }
490        assert (jRow==smallestIndex);
491      }
492#endif
493      //smallestIndex=pivotLinkedForwardsAddress_[-1];
494      updateColumnTransposeU ( regionSparse, smallestIndex
495#if ABC_SMALL<2
496                  ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
497#endif
498#if ABC_PARALLEL
499                    ,0
500#endif
501);
502    }
503  } else {
504#endif
505#if ABC_SMALL>=0
506    // No row copy check space
507    CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
508    CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
509    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
510    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
511    nInRow=replaceColumnU(regionSparse,deletedPosition,deletedColumns,pivotRow);
512#endif
513#if ABC_SMALL<2
514  }
515#endif
516  numberNonZero = regionSparse->getNumElements (  );
517  //printf("replace nels %d\n",numberNonZero);
518  CoinFactorizationDouble saveFromU = 0.0;
519
520  for (CoinBigIndex i = 0; i < numberInColumnU2; i++ ) {
521    CoinSimplexInt iRow = indexU2[i];
522    if ( iRow != pivotRow ) {
523      saveFromU -= elementU2[i] * region[iRow];
524    } else {
525      saveFromU += elementU2[i];
526    }
527  }       
528  //check accuracy
529  status = checkPivot(saveFromU,pivotCheck);
530  if (status==1&&!numberPivots_) {
531    printf("check status ok\n");
532    status=2;
533  }
534  pivotCheck=saveFromU;
535  return status;
536}
537//  replaceColumn.  Replaces one Column to basis
538//      returns 0=OK, 1=Probably OK, 2=singular, 3=no room
539//partial update already in U
540CoinSimplexInt
541CoinAbcTypeFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
542                                 CoinSimplexInt pivotRow,
543                                  CoinSimplexDouble pivotCheck ,
544                                  bool skipBtranU,
545                                   CoinSimplexDouble )
546{
547  assert (numberU_<=numberRowsExtra_);
548 
549#ifndef ALWAYS_SKIP_BTRAN
550  //return at once if too many iterations
551  if ( numberRowsExtra_ >= maximumRowsExtra_ ) {
552    return 5;
553  }       
554  if ( lengthAreaU_ < lastEntryByColumnU_ ) {
555    return 3;
556  }
557#endif 
558#ifndef ABC_USE_FUNCTION_POINTERS
559  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
560  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
561  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
562#endif
563  CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
564#if ABC_SMALL<2
565  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
566#endif
567#if ABC_SMALL<2
568  CoinSimplexInt * COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
569#endif
570  CoinSimplexInt *  COIN_RESTRICT permuteLookup = pivotColumnAddress_;
571  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
572 
573  //take out old pivot column
574 
575#ifndef ABC_USE_FUNCTION_POINTERS
576  totalElements_ -= numberInColumn[pivotRow];
577#else
578  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
579  extern scatterUpdate AbcScatterLowSubtract[9];
580  extern scatterUpdate AbcScatterHighSubtract[4];
581  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
582  CoinBigIndex startU = scatter[numberRows_].offset;
583  totalElements_ -= scatter[pivotRow].number;
584  CoinBigIndex saveStart = scatter[pivotRow].offset;
585  CoinBigIndex saveEnd = scatter[pivotRow].number;
586  scatter[pivotRow].offset=startU;
587  CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
588  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
589  CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
590#endif
591#ifndef ALWAYS_SKIP_BTRAN
592
593#ifndef ABC_USE_FUNCTION_POINTERS
594  CoinBigIndex saveStart = startColumnU[pivotRow];
595  CoinBigIndex saveEnd = saveStart + numberInColumn[pivotRow];
596#endif
597#endif
598  //get entries in row (pivot not stored)
599  CoinSimplexInt numberNonZero = 0;
600#if ABC_SMALL<2
601  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
602#endif
603#if CONVERTROW
604  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
605#if CONVERTROW>2
606  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
607#endif
608#endif
609#if ABC_SMALL<2
610  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
611  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
612#endif
613  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
614  CoinFactorizationDouble pivotValue = 1.0;
615  CoinSimplexInt status=0;
616#ifndef ABC_USE_FUNCTION_POINTERS
617  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
618#endif
619  CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
620
621#ifndef NDEBUG
622#if CONVERTROW
623#define CONVERTDEBUG
624#endif
625#endif
626  int nInRow;
627#if ABC_SMALL<2
628  CoinBigIndex start=0;
629  CoinBigIndex end=0;
630  if (gotUCopy()) {
631    start=startRowU[pivotRow];
632    nInRow=numberInRowU[pivotRow];
633    end= start + nInRow;
634  }
635#endif
636  CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
637#ifndef ALWAYS_SKIP_BTRAN
638  if (!skipBtranU) {
639  CoinFactorizationDouble oldPivot = pivotRegion[pivotRow];
640  // for accuracy check
641  pivotCheck = pivotCheck / oldPivot;
642#if ABC_SMALL<2
643  if (gotUCopy()) {
644    CoinSimplexInt smallestIndex=pivotLinkedForwardsAddress_[-1];
645    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
646    if (nInRow<10) {
647      CoinSimplexInt smallest=numberRowsExtra_;
648      for (CoinBigIndex i = start; i < end ; i ++ ) {
649        CoinSimplexInt jColumn = indexColumnU[i];
650        if (pivotRowForward[jColumn]<smallest) {
651          smallest=pivotRowForward[jColumn];
652          smallestIndex=jColumn;
653        }
654#ifndef ABC_USE_FUNCTION_POINTERS
655#ifdef CONVERTDEBUG
656        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
657        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
658#endif     
659        region[jColumn] = elementRowU[i];
660#else
661#ifdef CONVERTDEBUG
662        CoinBigIndex j = convertRowToColumn[i];
663        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
664        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
665#endif     
666        region[jColumn] = elementRowU[i];
667#endif
668        regionIndex[numberNonZero++] = jColumn;
669      }
670    } else {
671      for (CoinBigIndex i = start; i < end ; i ++ ) {
672        CoinSimplexInt jColumn = indexColumnU[i];
673#ifdef CONVERTDEBUG
674#ifdef ABC_USE_FUNCTION_POINTERS
675        CoinBigIndex j = convertRowToColumn[i];
676        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
677        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
678#else
679        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
680        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
681#endif
682#endif     
683        region[jColumn] = elementRowU[i];
684        regionIndex[numberNonZero++] = jColumn;
685      }
686    }
687    //do BTRAN - finding first one to use
688    regionSparse->setNumElements ( numberNonZero );
689    if (numberNonZero) {
690      assert (smallestIndex>=0);
691#ifndef NDEBUG
692      const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
693      CoinSimplexInt jRow=pivotLinked[-1];
694      if (jRow!=smallestIndex) {
695        while (jRow>=0) {
696          CoinFactorizationDouble pivotValue = region[jRow];
697          if (pivotValue)
698            break;
699          jRow=pivotLinked[jRow];
700        }
701        assert (jRow==smallestIndex);
702      }
703#endif
704      //smallestIndex=pivotLinkedForwardsAddress_[-1];
705      updateColumnTransposeU ( regionSparse, smallestIndex
706#if ABC_SMALL<2
707                  ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
708#endif
709#if ABC_PARALLEL
710                    ,0
711#endif
712);
713    }
714  } else {
715#endif
716#if ABC_SMALL>=0
717    // No row copy check space
718    if ( lengthR_ + numberRows_>= lengthAreaR_ ) {
719      //not enough room
720      return 3;
721    }       
722    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
723    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
724    nInRow=replaceColumnU(regionSparse,deletedPosition,deletedColumns,pivotRow);
725#endif
726#if ABC_SMALL<2
727  }
728#endif
729  }
730#endif
731  numberNonZero = regionSparse->getNumElements (  );
732  CoinFactorizationDouble saveFromU = 0.0;
733
734#ifndef ABC_USE_FUNCTION_POINTERS
735  CoinBigIndex startU = startColumnU[numberRows_];
736  startColumnU[pivotRow]=startU;
737  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
738  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
739  CoinSimplexInt numberInColumnU2 = numberInColumn[numberRows_];
740#else
741//CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
742#endif
743  instrument_do("CoinAbcFactorizationReplaceColumn",2*numberRows_+4*numberInColumnU2);
744#ifndef ALWAYS_SKIP_BTRAN
745  if (!skipBtranU) {
746  for (CoinBigIndex i = 0; i < numberInColumnU2; i++ ) {
747    CoinSimplexInt iRow = indexU2[i];
748    if ( iRow != pivotRow ) {
749      saveFromU -= elementU2[i] * region[iRow];
750    } else {
751      saveFromU += elementU2[i];
752    }
753  }       
754  //check accuracy
755  status = checkPivot(saveFromU,pivotCheck);
756  } else {
757    status=0;
758    saveFromU=pivotCheck;
759  }
760  //regionSparse->print();
761#endif
762  if (status>1||(status==1&&!numberPivots_)) {
763    // restore some things
764    //pivotRegion[pivotRow] = oldPivot;
765    CoinSimplexInt number = saveEnd-saveStart;
766    totalElements_ += number;
767    //numberInColumn[pivotRow]=number;
768    regionSparse->clear();
769    return status;
770  } else {
771    pivotValue = 1.0 / saveFromU;
772    // do what we would have done by now
773#if ABC_SMALL<2
774    if (gotUCopy()) {
775      for (CoinBigIndex i = start; i < end ; i ++ ) {
776        CoinSimplexInt jColumn = indexColumnU[i];
777#ifndef ABC_USE_FUNCTION_POINTERS
778        CoinBigIndex startColumn = startColumnU[jColumn];
779#if CONVERTROW
780        CoinBigIndex j = convertRowToColumn[i];
781#else
782        CoinBigIndex j;
783        int number = numberInColumn[jColumn];
784        for (j=0;j<number;j++) {
785          if (indexRowU[j+startColumn]==pivotRow)
786            break;
787        }
788        assert (j<number);
789        //assert (j==convertRowToColumn[i]); // temp
790#endif
791#ifdef CONVERTDEBUG
792        assert (fabs(elementU[j+startColumn]-elementRowU[i])<1.0e-4);
793#endif     
794#else
795        int number = scatter[jColumn].number;
796        CoinSimplexInt k = number-1;
797        scatter[jColumn].number=k;
798        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
799        CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+k+1);
800#if CONVERTROW
801        CoinSimplexInt j = convertRowToColumn[i];
802#else
803        CoinSimplexInt j;
804        for (j=0;j<number;j++) {
805          if (indices[j]==pivotRow)
806            break;
807        }
808        assert (j<number);
809        //assert (j==convertRowToColumn[i]); // temp
810#endif
811#ifdef CONVERTDEBUG
812        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
813#endif     
814#endif
815        // swap
816#ifndef ABC_USE_FUNCTION_POINTERS
817        CoinSimplexInt k = numberInColumn[jColumn]-1;
818        numberInColumn[jColumn]=k;
819        CoinBigIndex k2 = k+startColumn;
820        int kRow2=indexRowU[k2];
821        indexRowU[j+startColumn]=kRow2;
822        CoinFactorizationDouble value2=elementU[k2];
823        elementU[j+startColumn]=value2;
824#else
825        int kRow2=indices[k];
826        CoinFactorizationDouble value2=area[k];
827#if ABC_USE_FUNCTION_POINTERS
828        if (k<9) {
829          scatter[jColumn].functionPointer=AbcScatterLowSubtract[k];
830        } else {
831          scatter[jColumn].functionPointer=AbcScatterHighSubtract[k&3];
832        }
833#endif
834        // later more complicated swap to avoid move
835        indices[j]=kRow2;
836        area[j]=value2;
837        // move
838        indices-=2;
839        for (int i=0;i<k;i++)
840          indices[i]=indices[i+2];
841#endif
842        // move in row copy (slow - as other remark (but shouldn't need to))
843        CoinBigIndex start2 = startRowU[kRow2];
844        int n=numberInRowU[kRow2];
845        CoinBigIndex end2 = start2 + n;
846#ifndef NDEBUG
847        bool found=false;
848#endif
849        for (CoinBigIndex i2 = start2; i2 < end2 ; i2 ++ ) {
850          CoinSimplexInt iColumn2 = indexColumnU[i2];
851          if (jColumn==iColumn2) {
852#if CONVERTROW
853            convertRowToColumn[i2] = j;
854#endif
855#ifndef NDEBUG
856            found=true;
857#endif
858            break;
859          }
860        }
861#ifndef NDEBUG
862        assert(found);
863#endif
864      }
865    } else {
866#endif
867#if ABC_SMALL>=0
868      // no row copy
869      CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
870      CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
871      for (CoinSimplexInt i = 0; i < nInRow ; i ++ ) {
872        CoinBigIndex j = deletedPosition[i];
873        CoinSimplexInt jColumn = deletedColumns[i];
874        // swap
875        CoinSimplexInt k = numberInColumn[jColumn]-1;
876        numberInColumn[jColumn]=k;
877        CoinBigIndex k2 = k+startColumnU[jColumn];
878        int kRow2=indexRowU[k2];
879        indexRowU[j]=kRow2;
880        CoinFactorizationDouble value2=elementU[k2];
881        elementU[j]=value2;
882      }
883#endif
884#if ABC_SMALL<2
885    }
886#endif
887  }
888#if ABC_SMALL<2
889  if (gotUCopy()) {
890    // Now zero out column of U
891    //take out old pivot column
892#ifdef ABC_USE_FUNCTION_POINTERS
893    CoinFactorizationDouble * COIN_RESTRICT elementU = elementUColumnPlus+saveStart;
894    saveStart=0;
895    CoinSimplexInt * COIN_RESTRICT indexRowU = reinterpret_cast<CoinSimplexInt *>(elementU+saveEnd);
896#endif
897    for (CoinBigIndex i = saveStart; i < saveEnd ; i ++ ) {
898      //elementU[i] = 0.0;
899      // If too slow then reverse meaning of convertRowToColumn and use
900      int jRow=indexRowU[i];
901      CoinBigIndex start = startRowU[jRow];
902      CoinBigIndex end = start + numberInRowU[jRow];
903      for (CoinBigIndex j = start; j < end ; j ++ ) {
904        CoinSimplexInt jColumn = indexColumnU[j];
905        if (jColumn==pivotRow) {
906          // swap
907          numberInRowU[jRow]--;
908          elementRowU[j]=elementRowU[end-1];
909#if CONVERTROW
910          convertRowToColumn[j]=convertRowToColumn[end-1];
911#endif
912          indexColumnU[j]=indexColumnU[end-1];
913          break;
914        }
915      }
916    }   
917  }
918#endif   
919  //zero out pivot Row (before or after?)
920  //add to R
921  CoinBigIndex * COIN_RESTRICT startColumnR = startColumnRAddress_;
922  CoinBigIndex putR = lengthR_;
923  CoinSimplexInt number = numberR_;
924 
925  startColumnR[number] = putR;  //for luck and first time
926  number++;
927  //assert (startColumnR+number-firstCountAddress_<( CoinMax(5*numberRows_,2*numberRows_+2*maximumPivots_)+2));
928  startColumnR[number] = putR + numberNonZero;
929  numberR_ = number;
930  lengthR_ = putR + numberNonZero;
931  totalElements_ += numberNonZero;
932  if ( lengthR_ >= lengthAreaR_ ) {
933    //not enough room
934    regionSparse->clear();
935    return 3;
936  }       
937  if ( lengthU_+numberInColumnU2 >= lengthAreaU_ ) {
938    //not enough room
939    regionSparse->clear();
940    return 3;
941  }
942
943#if ABC_SMALL<2
944  CoinSimplexInt *  COIN_RESTRICT nextRow=NULL;
945  CoinSimplexInt *  COIN_RESTRICT lastRow=NULL;
946  CoinSimplexInt next;
947  CoinSimplexInt last;
948  if (gotUCopy()) {
949    //take out row
950    nextRow = nextRowAddress_;
951    lastRow = lastRowAddress_;
952    next = nextRow[pivotRow];
953    last = lastRow[pivotRow];
954   
955    nextRow[last] = next;
956    lastRow[next] = last;
957    numberInRowU[pivotRow]=0;
958  }
959#endif 
960  //put in pivot
961  //add row counts
962 
963  int n = 0;
964  for (CoinSimplexInt i = 0; i < numberInColumnU2; i++ ) {
965    CoinSimplexInt iRow = indexU2[i];
966    CoinFactorizationDouble value=elementU2[i];
967    assert (value);
968    if ( !TEST_LESS_THAN_TOLERANCE(value ) ) {
969      if ( iRow != pivotRow ) {
970        //modify by pivot
971        CoinFactorizationDouble value2 = value*pivotValue;
972        indexU2[n]=iRow;
973        elementU2[n++] = value2;
974#if ABC_SMALL<2
975        if (gotUCopy()) {
976          CoinSimplexInt next = nextRow[iRow];
977          CoinSimplexInt iNumberInRow = numberInRowU[iRow];
978          CoinBigIndex space;
979          CoinBigIndex put = startRowU[iRow] + iNumberInRow;
980          space = startRowU[next] - put;
981          if ( space <= 0 ) {
982            getRowSpaceIterate ( iRow, iNumberInRow + 4 );
983            put = startRowU[iRow] + iNumberInRow;
984          } else if (next==numberRows_) {
985            lastEntryByRowU_=put+1;
986          }     
987          indexColumnU[put] = pivotRow;
988#if CONVERTROW
989          convertRowToColumn[put] = i;
990#endif
991          elementRowU[put] = value2;
992          numberInRowU[iRow] = iNumberInRow + 1;
993        }
994#endif
995      } else {
996        //swap and save
997        indexU2[i]=indexU2[numberInColumnU2-1];
998        elementU2[i] = elementU2[numberInColumnU2-1];
999        numberInColumnU2--;
1000        i--;
1001      }
1002    } else {
1003      // swap
1004      indexU2[i]=indexU2[numberInColumnU2-1];
1005      elementU2[i] = elementU2[numberInColumnU2-1];
1006      numberInColumnU2--;
1007      i--;
1008    }       
1009  }       
1010  numberInColumnU2=n;
1011  //do permute
1012  permuteAddress_[numberRowsExtra_] = pivotRow;
1013  //and for safety
1014  permuteAddress_[numberRowsExtra_ + 1] = 0;
1015
1016  //pivotColumnAddress_[pivotRow] = numberRowsExtra_;
1017 
1018  numberU_++;
1019
1020#if ABC_SMALL<2
1021  if (gotUCopy()) {
1022    //in at end
1023    last = lastRow[numberRows_];
1024    nextRow[last] = pivotRow; //numberRowsExtra_;
1025    lastRow[numberRows_] = pivotRow; //numberRowsExtra_;
1026    lastRow[pivotRow] = last;
1027    nextRow[pivotRow] = numberRows_;
1028    startRowU[pivotRow] = lastEntryByRowU_;
1029  }
1030#endif
1031#ifndef ABC_USE_FUNCTION_POINTERS
1032  numberInColumn[pivotRow]=numberInColumnU2;
1033#else
1034#if ABC_USE_FUNCTION_POINTERS
1035  if (numberInColumnU2<9) {
1036    scatter[pivotRow].functionPointer=AbcScatterLowSubtract[numberInColumnU2];
1037  } else {
1038    scatter[pivotRow].functionPointer=AbcScatterHighSubtract[numberInColumnU2&3];
1039  }
1040#endif
1041  assert(scatter[pivotRow].offset==lastEntryByColumnUPlus_);
1042  //CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+lastEntryByColumnUPlus_;
1043  //CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+numberInColumnU2);
1044  if (numberInColumnU2<scatter[numberRows_].number) {
1045    int offset=2*(scatter[numberRows_].number-numberInColumnU2);
1046    indexU2 -= offset;
1047    //assert(indexU2-reinterpret_cast<int*>(area)>=0);
1048    for (int i=0;i<numberInColumnU2;i++)
1049      indexU2[i]=indexU2[i+offset];
1050  }
1051  scatter[pivotRow].number=numberInColumnU2;
1052  //CoinAbcMemcpy(indices,indexU2,numberInColumnU2);
1053  //CoinAbcMemcpy(area,elementU2,numberInColumnU2);
1054  lastEntryByColumnUPlus_ += (3*numberInColumnU2+1)>>1;
1055#endif
1056  totalElements_ += numberInColumnU2;
1057  lengthU_ += numberInColumnU2;
1058#if ABC_SMALL<2
1059  CoinSimplexInt * COIN_RESTRICT nextColumn = nextColumnAddress_;
1060  CoinSimplexInt * COIN_RESTRICT lastColumn = lastColumnAddress_;
1061  if (gotRCopy()) {
1062    //column in at beginning (as empty)
1063    next = nextColumn[maximumRowsExtra_];
1064    lastColumn[next] = numberRowsExtra_;
1065    nextColumn[maximumRowsExtra_] = numberRowsExtra_;
1066    nextColumn[numberRowsExtra_] = next;
1067    lastColumn[numberRowsExtra_] = maximumRowsExtra_;
1068  }
1069#endif
1070#if 0
1071  {
1072    int k=-1;
1073    for (int i=0;i<numberRows_;i++) {
1074      k=CoinMax(k,startRowU[i]+numberInRowU[i]);
1075    }
1076    assert (k==lastEntryByRowU_);
1077  }
1078#endif
1079  //pivotRegion[numberRowsExtra_] = pivotValue;
1080  pivotRegion[pivotRow] = pivotValue;
1081#ifndef ABC_USE_FUNCTION_POINTERS
1082  maximumU_ = CoinMax(maximumU_,startU+numberInColumnU2);
1083#else
1084  maximumU_ = CoinMax(maximumU_,lastEntryByColumnUPlus_);
1085#endif
1086  permuteLookup[pivotRow]=numberRowsExtra_;
1087  permuteLookup[numberRowsExtra_]=pivotRow;
1088  //pivotColumnAddress_[pivotRow]=numberRowsExtra_;
1089  //pivotColumnAddress_[numberRowsExtra_]=pivotRow;
1090  assert (pivotRow<numberRows_);
1091  // out of chain
1092  CoinSimplexInt iLast=pivotLinkedBackwardsAddress_[pivotRow];
1093  CoinSimplexInt iNext=pivotLinkedForwardsAddress_[pivotRow];
1094  assert (pivotRow==pivotLinkedForwardsAddress_[iLast]);
1095  assert (pivotRow==pivotLinkedBackwardsAddress_[iNext]);
1096  pivotLinkedForwardsAddress_[iLast]=iNext;
1097  pivotLinkedBackwardsAddress_[iNext]=iLast;
1098  if (pivotRow==lastSlack_) {
1099    lastSlack_ = iLast;
1100  }
1101  iLast=pivotLinkedBackwardsAddress_[numberRows_];
1102  assert (numberRows_==pivotLinkedForwardsAddress_[iLast]);
1103  pivotLinkedForwardsAddress_[iLast]=pivotRow;
1104  pivotLinkedBackwardsAddress_[pivotRow]=iLast;
1105  pivotLinkedBackwardsAddress_[numberRows_]=pivotRow;
1106  pivotLinkedForwardsAddress_[pivotRow]=numberRows_;
1107  assert (numberRows_>pivotLinkedForwardsAddress_[-1]);
1108  assert (pivotLinkedBackwardsAddress_[numberRows_]>=0);
1109  numberRowsExtra_++;
1110  numberGoodU_++;
1111  numberPivots_++;
1112  if ( numberRowsExtra_ > numberRows_ + 50 ) {
1113    CoinBigIndex extra = factorElements_ >> 1;
1114   
1115    if ( numberRowsExtra_ > numberRows_ + 100 + numberRows_ / 500 ) {
1116      if ( extra < 2 * numberRows_ ) {
1117        extra = 2 * numberRows_;
1118      }       
1119    } else {
1120      if ( extra < 5 * numberRows_ ) {
1121        extra = 5 * numberRows_;
1122      }       
1123    }       
1124    CoinBigIndex added = totalElements_ - factorElements_;
1125   
1126    if ( added > extra && added > ( factorElements_ ) << 1 && !status
1127         && 3*totalElements_ > 2*(lengthAreaU_+lengthAreaL_)) {
1128      status = 3;
1129      if ( messageLevel_ & 4 ) {
1130        std::cout << "Factorization has "<< totalElements_
1131                  << ", basis had " << factorElements_ <<std::endl;
1132      }
1133    }       
1134  }
1135#if ABC_SMALL<2
1136  if (gotRCopy()&&status<2) {
1137    //if (numberInColumnPlus&&status<2) {
1138    // we are going to put another copy of R in R
1139    CoinFactorizationDouble * COIN_RESTRICT elementRR = elementRAddress_ + lengthAreaR_;
1140    CoinSimplexInt * COIN_RESTRICT indexRowRR = indexRowRAddress_ + lengthAreaR_;
1141    CoinBigIndex * COIN_RESTRICT startRR = startColumnRAddress_+maximumPivots_+1;
1142    CoinSimplexInt pivotRow = numberRowsExtra_-1;
1143    for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
1144      CoinSimplexInt jRow = regionIndex[i];
1145      CoinFactorizationDouble value=region[jRow];
1146      elementR[putR] = value;
1147      indexRowR[putR] = jRow;
1148      putR++;
1149      region[jRow]=0.0;
1150      int iRow = permuteLookup[jRow];
1151      next = nextColumn[iRow];
1152      CoinBigIndex space;
1153      if (next!=maximumRowsExtra_)
1154        space = startRR[next]-startRR[iRow];
1155      else
1156        space = lengthAreaR_-startRR[iRow];
1157      CoinSimplexInt numberInR = numberInColumnPlus[iRow];
1158      if (space>numberInR) {
1159        // there is space
1160        CoinBigIndex  put=startRR[iRow]+numberInR;
1161        numberInColumnPlus[iRow]=numberInR+1;
1162        indexRowRR[put]=pivotRow;
1163        elementRR[put]=value;
1164        //add 4 for luck
1165        if (next==maximumRowsExtra_)
1166          startRR[maximumRowsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaR_);
1167      } else {
1168        // no space - do we shuffle?
1169        if (!getColumnSpaceIterateR(iRow,value,pivotRow)) {
1170          // printf("Need more space for R\n");
1171          numberInColumnPlus_.conditionalDelete();
1172          numberInColumnPlusAddress_=NULL;
1173          setNoGotRCopy();
1174          regionSparse->clear();
1175#if ABC_SMALL<0
1176          status=3;
1177#endif
1178          break;
1179        }
1180      }
1181    }       
1182  } else {
1183#endif
1184    for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
1185      CoinSimplexInt iRow = regionIndex[i];
1186      elementR[putR] = region[iRow];
1187      region[iRow]=0.0;
1188      indexRowR[putR] = iRow;
1189      putR++;
1190    }       
1191#if ABC_SMALL<2
1192  }
1193#endif
1194#ifdef ABC_USE_FUNCTION_POINTERS
1195#if 0
1196  {
1197    CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1198    CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
1199    CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
1200    CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1201    CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
1202    scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
1203    extern scatterUpdate AbcScatterLowSubtract[9];
1204    extern scatterUpdate AbcScatterHighSubtract[4];
1205    printf("=============================================start\n");
1206    for (int iRow=0;iRow<numberRows_;iRow++) {
1207      int number=scatter[iRow].number;
1208      CoinBigIndex offset = scatter[iRow].offset;
1209      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+offset;
1210      CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+number);
1211      CoinSimplexInt * COIN_RESTRICT indices2 = indexRow+startColumnU[iRow];
1212      CoinFactorizationDouble * COIN_RESTRICT area2 = element+startColumnU[iRow];
1213      //assert (number==numberInColumn[iRow]);
1214#if ABC_USE_FUNCTION_POINTERS
1215      if (number<9)
1216        assert (scatter[iRow].functionPointer==AbcScatterLowSubtract[number]);
1217      else
1218        assert (scatter[iRow].functionPointer==AbcScatterHighSubtract[number&3]);
1219#endif
1220      if (number)
1221        printf("Row %d %d els\n",iRow,number);
1222      int n=0;
1223      for (int i=0;i<number;i++) {
1224        printf("(%d,%g) ",indices[i],area[i]);
1225        n++;
1226        if (n==10) {
1227          n=0;
1228          printf("\n");
1229        }
1230      }
1231      if (n)
1232        printf("\n");
1233#if 0
1234      for (int i=0;i<number;i++) {
1235        assert (indices[i]==indices2[i]);
1236        assert (area[i]==area2[i]);
1237      }
1238#endif
1239    }
1240    printf("=============================================end\n");
1241  }
1242#endif
1243#endif
1244  regionSparse->setNumElements(0);
1245  //for (int i=0;i<numberRows_;i++)
1246  //assert(permuteLookupAddress_[i]==pivotColumnAddress_[i]);
1247  return status;
1248}
1249#endif
1250/* Checks if can replace one Column to basis,
1251   returns update alpha
1252   Fills in region for use later
1253   partial update already in U */
1254#ifdef ABC_LONG_FACTORIZATION
1255  long
1256#endif
1257double
1258CoinAbcTypeFactorization::checkReplacePart1 (CoinIndexedVector * regionSparse,
1259                                              int pivotRow)
1260{
1261#ifdef ABC_USE_FUNCTION_POINTERS
1262  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
1263  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1264  CoinBigIndex startU = scatter[numberRows_].offset;
1265  CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
1266  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
1267  //CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
1268#else
1269  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
1270  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1271  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
1272  CoinBigIndex startU = startColumnU[numberRows_];
1273  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
1274  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
1275  CoinSimplexInt numberInColumnU2 = numberInColumn[numberRows_];
1276#endif
1277#if ABC_SMALL<2
1278  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
1279#endif
1280  //zeroed out region
1281  toLongArray(regionSparse,3);
1282  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
1283  //CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1284
1285#ifndef ABC_USE_FUNCTION_POINTERS
1286  CoinBigIndex saveStart = startColumnU[pivotRow];
1287  CoinBigIndex saveEnd = saveStart + numberInColumn[pivotRow];
1288#endif
1289  //get entries in row (pivot not stored)
1290  CoinSimplexInt numberNonZero = 0;
1291#if ABC_SMALL<2
1292  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
1293#endif
1294#if CONVERTROW
1295#ifndef NDEBUG
1296#define CONVERTDEBUG
1297#endif
1298#ifdef CONVERTDEBUG
1299  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
1300#endif
1301#if CONVERTROW>2
1302  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
1303#endif
1304#endif
1305#if ABC_SMALL<2
1306  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
1307  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
1308#endif
1309  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1310  //CoinSimplexInt status=0;
1311#ifndef ABC_USE_FUNCTION_POINTERS
1312  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
1313#endif
1314
1315  int nInRow;
1316#if ABC_SMALL<2
1317  CoinBigIndex start=0;
1318  CoinBigIndex end=0;
1319  if (gotUCopy()) {
1320    start=startRowU[pivotRow];
1321    nInRow=numberInRowU[pivotRow];
1322    end= start + nInRow;
1323    CoinSimplexInt smallestIndex=pivotLinkedForwardsAddress_[-1];
1324    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
1325    if (nInRow<10) {
1326      CoinSimplexInt smallest=numberRowsExtra_;
1327      for (CoinBigIndex i = start; i < end ; i ++ ) {
1328        CoinSimplexInt jColumn = indexColumnU[i];
1329        if (pivotRowForward[jColumn]<smallest) {
1330          smallest=pivotRowForward[jColumn];
1331          smallestIndex=jColumn;
1332        }
1333#ifndef ABC_USE_FUNCTION_POINTERS
1334#ifdef CONVERTDEBUG
1335        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
1336        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
1337#endif     
1338        region[jColumn] = elementRowU[i];
1339#else
1340#ifdef CONVERTDEBUG
1341        CoinBigIndex j = convertRowToColumn[i];
1342        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1343        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
1344#endif     
1345        region[jColumn] = elementRowU[i];
1346#endif
1347        regionIndex[numberNonZero++] = jColumn;
1348      }
1349    } else {
1350      for (CoinBigIndex i = start; i < end ; i ++ ) {
1351        CoinSimplexInt jColumn = indexColumnU[i];
1352#ifdef CONVERTDEBUG
1353#ifdef ABC_USE_FUNCTION_POINTERS
1354        CoinBigIndex j = convertRowToColumn[i];
1355        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1356        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
1357#else
1358        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
1359        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
1360#endif
1361#endif     
1362        region[jColumn] = elementRowU[i];
1363        regionIndex[numberNonZero++] = jColumn;
1364      }
1365    }
1366    //do BTRAN - finding first one to use
1367    regionSparse->setNumElements ( numberNonZero );
1368    if (numberNonZero) {
1369      assert (smallestIndex>=0);
1370#ifndef NDEBUG
1371      const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
1372      CoinSimplexInt jRow=pivotLinked[-1];
1373      if (jRow!=smallestIndex) {
1374        while (jRow>=0) {
1375          CoinFactorizationDouble pivotValue = region[jRow];
1376          if (pivotValue)
1377            break;
1378          jRow=pivotLinked[jRow];
1379        }
1380        assert (jRow==smallestIndex);
1381      }
1382#endif
1383#if ABC_PARALLEL==0
1384      //smallestIndex=pivotLinkedForwardsAddress_[-1];
1385      updateColumnTransposeU ( regionSparse, smallestIndex
1386#if ABC_SMALL<2
1387              ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
1388#endif
1389                               );
1390#else
1391      assert (FACTOR_CPU>3);
1392      updateColumnTransposeU ( regionSparse, smallestIndex
1393#if ABC_SMALL<2
1394              ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
1395#endif
1396                               ,3
1397                               );
1398#endif
1399    }
1400  } else {
1401#endif
1402#if ABC_SMALL>=0
1403    // No row copy check space
1404    CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
1405    CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
1406    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
1407    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
1408    nInRow=replaceColumnU(regionSparse,deletedPosition,deletedColumns,pivotRow);
1409#endif
1410#if ABC_SMALL<2
1411  }
1412#endif
1413  if ( lengthR_+numberRows_ >= lengthAreaR_ ) {
1414    //not enough room
1415    return 0.0;
1416  }       
1417#ifdef ABC_USE_FUNCTION_POINTERS
1418  CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
1419#endif
1420  if ( lengthU_+numberInColumnU2 >= lengthAreaU_ ) {
1421    //not enough room
1422    return 0.0;
1423  }
1424  //CoinSimplexDouble * COIN_RESTRICT region = regionSparse->denseVector (  );
1425  //CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1426  //CoinSimplexInt status=0;
1427
1428  //int numberNonZero = regionSparse->getNumElements (  );
1429  //printf("replace nels %d\n",numberNonZero);
1430  CoinFactorizationDouble saveFromU = 0.0;
1431
1432  for (CoinBigIndex i = 0; i < numberInColumnU2; i++ ) {
1433    CoinSimplexInt iRow = indexU2[i];
1434    if ( iRow != pivotRow ) {
1435      saveFromU -= elementU2[i] * region[iRow];
1436    } else {
1437      saveFromU += elementU2[i];
1438    }
1439  }       
1440  return saveFromU;
1441}
1442/* Checks if can replace one Column to basis,
1443   returns update alpha
1444   Fills in region for use later
1445   partial update in vector */
1446#ifdef ABC_LONG_FACTORIZATION
1447  long
1448#endif
1449double
1450CoinAbcTypeFactorization::checkReplacePart1 (CoinIndexedVector * regionSparse,
1451                                             CoinIndexedVector * partialUpdate,
1452                                              int pivotRow)
1453{
1454#ifdef ABC_USE_FUNCTION_POINTERS
1455  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
1456  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1457  CoinBigIndex startU = scatter[numberRows_].offset;
1458#else
1459  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
1460  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1461  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
1462#endif
1463  CoinSimplexInt * COIN_RESTRICT indexU2 = partialUpdate->getIndices();
1464  CoinFactorizationDouble * COIN_RESTRICT elementU2 = denseVector(partialUpdate);
1465  CoinSimplexInt numberInColumnU2 = partialUpdate->getNumElements();
1466#if ABC_SMALL<2
1467  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
1468#endif
1469  //zeroed out region
1470  toLongArray(regionSparse,3);
1471  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
1472  //CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1473
1474#ifndef ABC_USE_FUNCTION_POINTERS
1475  CoinBigIndex saveStart = startColumnU[pivotRow];
1476  CoinBigIndex saveEnd = saveStart + numberInColumn[pivotRow];
1477#endif
1478  //get entries in row (pivot not stored)
1479  CoinSimplexInt numberNonZero = 0;
1480#if ABC_SMALL<2
1481  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
1482#endif
1483#if CONVERTROW
1484#ifndef NDEBUG
1485#define CONVERTDEBUG
1486#endif
1487#ifdef CONVERTDEBUG
1488  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
1489#endif
1490#if CONVERTROW>2
1491  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
1492#endif
1493#endif
1494#if ABC_SMALL<2
1495  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
1496  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
1497#endif
1498  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1499  //CoinSimplexInt status=0;
1500#ifndef ABC_USE_FUNCTION_POINTERS
1501  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
1502#endif
1503
1504  int nInRow;
1505#if ABC_SMALL<2
1506  CoinBigIndex start=0;
1507  CoinBigIndex end=0;
1508  if (gotUCopy()) {
1509    start=startRowU[pivotRow];
1510    nInRow=numberInRowU[pivotRow];
1511    end= start + nInRow;
1512    CoinSimplexInt smallestIndex=pivotLinkedForwardsAddress_[-1];
1513    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
1514    if (nInRow<10) {
1515      CoinSimplexInt smallest=numberRowsExtra_;
1516      for (CoinBigIndex i = start; i < end ; i ++ ) {
1517        CoinSimplexInt jColumn = indexColumnU[i];
1518        if (pivotRowForward[jColumn]<smallest) {
1519          smallest=pivotRowForward[jColumn];
1520          smallestIndex=jColumn;
1521        }
1522#ifndef ABC_USE_FUNCTION_POINTERS
1523#ifdef CONVERTDEBUG
1524        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
1525        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
1526#endif     
1527        region[jColumn] = elementRowU[i];
1528#else
1529#ifdef CONVERTDEBUG
1530        CoinBigIndex j = convertRowToColumn[i];
1531        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1532        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
1533#endif     
1534        region[jColumn] = elementRowU[i];
1535#endif
1536        regionIndex[numberNonZero++] = jColumn;
1537      }
1538    } else {
1539      for (CoinBigIndex i = start; i < end ; i ++ ) {
1540        CoinSimplexInt jColumn = indexColumnU[i];
1541#ifdef CONVERTDEBUG
1542#ifdef ABC_USE_FUNCTION_POINTERS
1543        CoinBigIndex j = convertRowToColumn[i];
1544        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1545        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
1546#else
1547        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
1548        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
1549#endif
1550#endif     
1551        region[jColumn] = elementRowU[i];
1552        regionIndex[numberNonZero++] = jColumn;
1553      }
1554    }
1555    //do BTRAN - finding first one to use
1556    regionSparse->setNumElements ( numberNonZero );
1557    if (numberNonZero) {
1558      assert (smallestIndex>=0);
1559#ifndef NDEBUG
1560      const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
1561      CoinSimplexInt jRow=pivotLinked[-1];
1562      if (jRow!=smallestIndex) {
1563        while (jRow>=0) {
1564          CoinFactorizationDouble pivotValue = region[jRow];
1565          if (pivotValue)
1566            break;
1567          jRow=pivotLinked[jRow];
1568        }
1569        assert (jRow==smallestIndex);
1570      }
1571#endif
1572#if ABC_PARALLEL==0
1573      //smallestIndex=pivotLinkedForwardsAddress_[-1];
1574      updateColumnTransposeU ( regionSparse, smallestIndex
1575#if ABC_SMALL<2
1576              ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
1577#endif
1578                               );
1579#else
1580      assert (FACTOR_CPU>3);
1581      updateColumnTransposeU ( regionSparse, smallestIndex
1582#if ABC_SMALL<2
1583              ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
1584#endif
1585                               ,3
1586                               );
1587#endif
1588    }
1589  } else {
1590#endif
1591#if ABC_SMALL>=0
1592    // No row copy check space
1593    CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
1594    CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
1595    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
1596    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
1597    nInRow=replaceColumnU(regionSparse,deletedPosition,deletedColumns,pivotRow);
1598#endif
1599#if ABC_SMALL<2
1600  }
1601#endif
1602  if ( lengthR_+numberRows_ >= lengthAreaR_ ) {
1603    //not enough room
1604    return 0.0;
1605  }       
1606  if ( lengthU_+numberInColumnU2 >= lengthAreaU_ ) {
1607    //not enough room
1608    partialUpdate->clear();
1609    return 0.0;
1610  }
1611  //CoinSimplexDouble * COIN_RESTRICT region = regionSparse->denseVector (  );
1612  //CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1613  //CoinSimplexInt status=0;
1614
1615  //int numberNonZero = regionSparse->getNumElements (  );
1616  //printf("replace nels %d\n",numberNonZero);
1617  CoinFactorizationDouble saveFromU = 0.0;
1618
1619  for (CoinBigIndex i = 0; i < numberInColumnU2; i++ ) {
1620    CoinSimplexInt iRow = indexU2[i];
1621    if ( iRow != pivotRow ) {
1622      saveFromU -= elementU2[iRow] * region[iRow];
1623    } else {
1624      saveFromU += elementU2[iRow];
1625    }
1626  }       
1627  return saveFromU;
1628}
1629#ifdef MOVE_REPLACE_PART1A
1630/* Checks if can replace one Column in basis,
1631   returns update alpha
1632   Fills in region for use later
1633   partial update already in U */
1634void
1635CoinAbcTypeFactorization::checkReplacePart1a (CoinIndexedVector * regionSparse,
1636                                              int pivotRow)
1637{
1638#ifdef ABC_USE_FUNCTION_POINTERS
1639  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
1640  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1641  CoinBigIndex startU = scatter[numberRows_].offset;
1642  CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
1643  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
1644  //CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
1645#else
1646  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
1647  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1648  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
1649  CoinBigIndex startU = startColumnU[numberRows_];
1650  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
1651  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
1652  CoinSimplexInt numberInColumnU2 = numberInColumn[numberRows_];
1653#endif
1654#if ABC_SMALL<2
1655  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
1656#endif
1657  //zeroed out region
1658  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
1659  //CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1660
1661#ifndef ABC_USE_FUNCTION_POINTERS
1662  CoinBigIndex saveStart = startColumnU[pivotRow];
1663  CoinBigIndex saveEnd = saveStart + numberInColumn[pivotRow];
1664#endif
1665  //get entries in row (pivot not stored)
1666  CoinSimplexInt numberNonZero = 0;
1667#if ABC_SMALL<2
1668  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
1669#endif
1670#if CONVERTROW
1671#ifndef NDEBUG
1672#define CONVERTDEBUG
1673#endif
1674#ifdef CONVERTDEBUG
1675  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
1676#endif
1677#if CONVERTROW>2
1678  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
1679#endif
1680#endif
1681#if ABC_SMALL<2
1682  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
1683  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
1684#endif
1685  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1686  //CoinSimplexInt status=0;
1687#ifndef ABC_USE_FUNCTION_POINTERS
1688  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
1689#endif
1690
1691  int nInRow;
1692#if ABC_SMALL<2
1693  CoinBigIndex start=0;
1694  CoinBigIndex end=0;
1695  if (gotUCopy()) {
1696    start=startRowU[pivotRow];
1697    nInRow=numberInRowU[pivotRow];
1698    end= start + nInRow;
1699    CoinSimplexInt smallestIndex=pivotLinkedForwardsAddress_[-1];
1700    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
1701    if (nInRow<10) {
1702      CoinSimplexInt smallest=numberRowsExtra_;
1703      for (CoinBigIndex i = start; i < end ; i ++ ) {
1704        CoinSimplexInt jColumn = indexColumnU[i];
1705        if (pivotRowForward[jColumn]<smallest) {
1706          smallest=pivotRowForward[jColumn];
1707          smallestIndex=jColumn;
1708        }
1709#ifndef ABC_USE_FUNCTION_POINTERS
1710#ifdef CONVERTDEBUG
1711        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
1712        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
1713#endif     
1714        region[jColumn] = elementRowU[i];
1715#else
1716#ifdef CONVERTDEBUG
1717        CoinBigIndex j = convertRowToColumn[i];
1718        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1719        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
1720#endif     
1721        region[jColumn] = elementRowU[i];
1722#endif
1723        regionIndex[numberNonZero++] = jColumn;
1724      }
1725    } else {
1726      for (CoinBigIndex i = start; i < end ; i ++ ) {
1727        CoinSimplexInt jColumn = indexColumnU[i];
1728#ifdef CONVERTDEBUG
1729#ifdef ABC_USE_FUNCTION_POINTERS
1730        CoinBigIndex j = convertRowToColumn[i];
1731        CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1732        assert (fabs(area[j]-elementRowU[i])<1.0e-4);
1733#else
1734        CoinBigIndex j = convertRowToColumn[i]+startColumnU[jColumn];
1735        assert (fabs(elementU[j]-elementRowU[i])<1.0e-4);
1736#endif
1737#endif     
1738        region[jColumn] = elementRowU[i];
1739        regionIndex[numberNonZero++] = jColumn;
1740      }
1741    }
1742    //do BTRAN - finding first one to use
1743    regionSparse->setNumElements ( numberNonZero );
1744    if (numberNonZero) {
1745      assert (smallestIndex>=0);
1746#ifndef NDEBUG
1747      const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
1748      CoinSimplexInt jRow=pivotLinked[-1];
1749      if (jRow!=smallestIndex) {
1750        while (jRow>=0) {
1751          CoinFactorizationDouble pivotValue = region[jRow];
1752          if (pivotValue)
1753            break;
1754          jRow=pivotLinked[jRow];
1755        }
1756        assert (jRow==smallestIndex);
1757      }
1758#endif
1759#if ABC_PARALLEL==0
1760      //smallestIndex=pivotLinkedForwardsAddress_[-1];
1761      updateColumnTransposeU ( regionSparse, smallestIndex
1762#if ABC_SMALL<2
1763              ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
1764#endif
1765                               );
1766#else
1767      assert (FACTOR_CPU>3);
1768      updateColumnTransposeU ( regionSparse, smallestIndex
1769#if ABC_SMALL<2
1770              ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
1771#endif
1772                               ,3
1773                               );
1774#endif
1775    }
1776  } else {
1777#endif
1778#if ABC_SMALL>=0
1779    // No row copy check space
1780    CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
1781    CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
1782    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
1783    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
1784    nInRow=replaceColumnU(regionSparse,deletedPosition,deletedColumns,pivotRow);
1785#endif
1786#if ABC_SMALL<2
1787  }
1788#endif
1789}
1790#ifdef ABC_LONG_FACTORIZATION
1791  long
1792#endif
1793double 
1794CoinAbcTypeFactorization::checkReplacePart1b (CoinIndexedVector * regionSparse,
1795                                              int pivotRow)
1796{
1797  if ( lengthR_+numberRows_ >= lengthAreaR_ ) {
1798    //not enough room
1799    return 0.0;
1800  }       
1801#ifdef ABC_USE_FUNCTION_POINTERS
1802  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
1803  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1804  CoinBigIndex startU = scatter[numberRows_].offset;
1805  CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
1806  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
1807  CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
1808#else
1809  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
1810  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1811  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
1812  CoinBigIndex startU = startColumnU[numberRows_];
1813  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
1814  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
1815  CoinSimplexInt numberInColumnU2 = numberInColumn[numberRows_];
1816#endif
1817  if ( lengthU_+numberInColumnU2 >= lengthAreaU_ ) {
1818    //not enough room
1819    return 0.0;
1820  }
1821  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
1822  //CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1823  //CoinSimplexInt status=0;
1824
1825  //int numberNonZero = regionSparse->getNumElements (  );
1826  //printf("replace nels %d\n",numberNonZero);
1827  CoinFactorizationDouble saveFromU = 0.0;
1828
1829  for (CoinBigIndex i = 0; i < numberInColumnU2; i++ ) {
1830    CoinSimplexInt iRow = indexU2[i];
1831    if ( iRow != pivotRow ) {
1832      saveFromU -= elementU2[i] * region[iRow];
1833    } else {
1834      saveFromU += elementU2[i];
1835    }
1836  }       
1837  return saveFromU;
1838}
1839#endif
1840#include "AbcSimplex.hpp"
1841/* Checks if can replace one Column to basis,
1842   returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots */
1843int 
1844CoinAbcTypeFactorization::checkReplacePart2 ( int pivotRow,
1845                                              CoinSimplexDouble /*btranAlpha*/, CoinSimplexDouble ftranAlpha, 
1846#ifdef ABC_LONG_FACTORIZATION
1847                                              long
1848#endif
1849                                              double ftAlpha,
1850                                              double /*acceptablePivot*/)
1851{
1852  if ( lengthR_+numberRows_ >= lengthAreaR_ ) {
1853    //not enough room
1854    return 3;
1855  }       
1856  CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1857  CoinFactorizationDouble oldPivot = pivotRegion[pivotRow];
1858  // for accuracy check
1859  CoinFactorizationDouble pivotCheck = ftranAlpha / oldPivot;
1860  //check accuracy
1861  int status = checkPivot(ftAlpha,pivotCheck);
1862  if (status==1&&!numberPivots_) {
1863    printf("check status ok\n");
1864    status=2;
1865  }
1866  return status;
1867}
1868/* Replaces one Column to basis,
1869   partial update already in U */
1870void 
1871CoinAbcTypeFactorization::replaceColumnPart3 ( const AbcSimplex * /*model*/,
1872                                          CoinIndexedVector * regionSparse,
1873                                               CoinIndexedVector * /*tableauColumn*/,
1874                                          int pivotRow,
1875#ifdef ABC_LONG_FACTORIZATION
1876                                               long
1877#endif
1878                                          double alpha )
1879{
1880  assert (numberU_<=numberRowsExtra_);
1881#ifndef ABC_USE_FUNCTION_POINTERS
1882  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
1883  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1884  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
1885#endif
1886  CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
1887#if ABC_SMALL<2
1888  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
1889#endif
1890#if ABC_SMALL<2
1891  CoinSimplexInt * COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
1892#endif
1893  CoinSimplexInt *  COIN_RESTRICT permuteLookup = pivotColumnAddress_;
1894  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
1895 
1896  //take out old pivot column
1897 
1898#ifndef ABC_USE_FUNCTION_POINTERS
1899  totalElements_ -= numberInColumn[pivotRow];
1900#else
1901  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
1902#if ABC_USE_FUNCTION_POINTERS
1903  extern scatterUpdate AbcScatterLowSubtract[9];
1904  extern scatterUpdate AbcScatterHighSubtract[4];
1905#endif
1906  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1907  CoinBigIndex startU = scatter[numberRows_].offset;
1908  totalElements_ -= scatter[pivotRow].number;
1909  CoinBigIndex saveStart = scatter[pivotRow].offset;
1910  CoinBigIndex saveEnd = scatter[pivotRow].number;
1911  scatter[pivotRow].offset=startU;
1912  CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
1913  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
1914  CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
1915#endif
1916  //get entries in row (pivot not stored)
1917  CoinSimplexInt numberNonZero = 0;
1918#if ABC_SMALL<2
1919  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
1920#endif
1921#if CONVERTROW
1922  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
1923#if CONVERTROW>2
1924  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
1925#endif
1926#endif
1927#if ABC_SMALL<2
1928  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
1929  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
1930#endif
1931  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
1932  CoinFactorizationDouble pivotValue = 1.0;
1933  CoinSimplexInt status=0;
1934#ifndef ABC_USE_FUNCTION_POINTERS
1935  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
1936#endif
1937  CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
1938
1939#ifndef NDEBUG
1940#if CONVERTROW
1941#define CONVERTDEBUG
1942#endif
1943#endif
1944  int nInRow=9999999; // ? what if ABC_SMALL==0 or ABC_SMALL==1
1945#if ABC_SMALL<2
1946  CoinBigIndex start=0;
1947  CoinBigIndex end=0;
1948  if (gotUCopy()) {
1949    start=startRowU[pivotRow];
1950    nInRow=numberInRowU[pivotRow];
1951    end= start + nInRow;
1952  }
1953#endif
1954  CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1955  numberNonZero = regionSparse->getNumElements (  );
1956  //CoinFactorizationDouble saveFromU = 0.0;
1957
1958#ifndef ABC_USE_FUNCTION_POINTERS
1959  CoinBigIndex startU = startColumnU[numberRows_];
1960  startColumnU[pivotRow]=startU;
1961  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
1962  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
1963  CoinSimplexInt numberInColumnU2 = numberInColumn[numberRows_];
1964#else
1965//CoinSimplexInt numberInColumnU2 = scatter[numberRows_].number;
1966#endif
1967  instrument_do("CoinAbcFactorizationReplaceColumn",2*numberRows_+4*numberInColumnU2);
1968  pivotValue = 1.0 / alpha;
1969#if ABC_SMALL<2
1970  if (gotUCopy()) {
1971    for (CoinBigIndex i = start; i < end ; i ++ ) {
1972      CoinSimplexInt jColumn = indexColumnU[i];
1973#ifndef ABC_USE_FUNCTION_POINTERS
1974      CoinBigIndex startColumn = startColumnU[jColumn];
1975#if CONVERTROW
1976      CoinBigIndex j = convertRowToColumn[i];
1977#else
1978      CoinBigIndex j;
1979      int number = numberInColumn[jColumn];
1980      for (j=0;j<number;j++) {
1981        if (indexRowU[j+startColumn]==pivotRow)
1982          break;
1983      }
1984      assert (j<number);
1985      //assert (j==convertRowToColumn[i]); // temp
1986#endif
1987#ifdef CONVERTDEBUG
1988      assert (fabs(elementU[j+startColumn]-elementRowU[i])<1.0e-4);
1989#endif     
1990#else
1991      int number = scatter[jColumn].number;
1992      CoinSimplexInt k = number-1;
1993      scatter[jColumn].number=k;
1994      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
1995      CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+k+1);
1996#if CONVERTROW
1997      CoinSimplexInt j = convertRowToColumn[i];
1998#else
1999      CoinSimplexInt j;
2000      for (j=0;j<number;j++) {
2001        if (indices[j]==pivotRow)
2002          break;
2003      }
2004      assert (j<number);
2005      //assert (j==convertRowToColumn[i]); // temp
2006#endif
2007#ifdef CONVERTDEBUG
2008      assert (fabs(area[j]-elementRowU[i])<1.0e-4);
2009#endif     
2010#endif
2011      // swap
2012#ifndef ABC_USE_FUNCTION_POINTERS
2013      CoinSimplexInt k = numberInColumn[jColumn]-1;
2014      numberInColumn[jColumn]=k;
2015      CoinBigIndex k2 = k+startColumn;
2016      int kRow2=indexRowU[k2];
2017      indexRowU[j+startColumn]=kRow2;
2018      CoinFactorizationDouble value2=elementU[k2];
2019      elementU[j+startColumn]=value2;
2020#else
2021      int kRow2=indices[k];
2022      CoinFactorizationDouble value2=area[k];
2023#if ABC_USE_FUNCTION_POINTERS
2024      if (k<9) {
2025        scatter[jColumn].functionPointer=AbcScatterLowSubtract[k];
2026      } else {
2027        scatter[jColumn].functionPointer=AbcScatterHighSubtract[k&3];
2028      }
2029#endif
2030      // later more complicated swap to avoid move (don't think we can!)
2031      indices[j]=kRow2;
2032      area[j]=value2;
2033      // move
2034      indices-=2;
2035      for (int i=0;i<k;i++)
2036        indices[i]=indices[i+2];
2037#endif
2038      // move in row copy (slow - as other remark (but shouldn't need to))
2039      CoinBigIndex start2 = startRowU[kRow2];
2040      int n=numberInRowU[kRow2];
2041      CoinBigIndex end2 = start2 + n;
2042#ifndef NDEBUG
2043      bool found=false;
2044#endif
2045      for (CoinBigIndex i2 = start2; i2 < end2 ; i2 ++ ) {
2046        CoinSimplexInt iColumn2 = indexColumnU[i2];
2047        if (jColumn==iColumn2) {
2048#if CONVERTROW
2049          convertRowToColumn[i2] = j;
2050#endif
2051#ifndef NDEBUG
2052          found=true;
2053#endif
2054          break;
2055        }
2056      }
2057#ifndef NDEBUG
2058      assert(found);
2059#endif
2060    }
2061  } else {
2062#endif
2063#if ABC_SMALL>=0
2064    assert(nInRow!=9999999); // ? what if ABC_SMALL==0 or ABC_SMALL==1
2065    // no row copy
2066    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
2067    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
2068    for (CoinSimplexInt i = 0; i < nInRow ; i ++ ) {
2069      CoinBigIndex j = deletedPosition[i];
2070      CoinSimplexInt jColumn = deletedColumns[i];
2071      // swap
2072      CoinSimplexInt k = numberInColumn[jColumn]-1;
2073      numberInColumn[jColumn]=k;
2074      CoinBigIndex k2 = k+startColumnU[jColumn];
2075      int kRow2=indexRowU[k2];
2076      indexRowU[j]=kRow2;
2077      CoinFactorizationDouble value2=elementU[k2];
2078      elementU[j]=value2;
2079    }
2080#endif
2081#if ABC_SMALL<2
2082  }
2083#endif
2084 #if ABC_SMALL<2
2085  if (gotUCopy()) {
2086    // Now zero out column of U
2087    //take out old pivot column
2088#ifdef ABC_USE_FUNCTION_POINTERS
2089    CoinFactorizationDouble * COIN_RESTRICT elementU = elementUColumnPlus+saveStart;
2090    saveStart=0;
2091    CoinSimplexInt * COIN_RESTRICT indexRowU = reinterpret_cast<CoinSimplexInt *>(elementU+saveEnd);
2092#endif
2093    for (CoinBigIndex i = saveStart; i < saveEnd ; i ++ ) {
2094      //elementU[i] = 0.0;
2095      // If too slow then reverse meaning of convertRowToColumn and use
2096      int jRow=indexRowU[i];
2097      CoinBigIndex start = startRowU[jRow];
2098      CoinBigIndex end = start + numberInRowU[jRow];
2099      for (CoinBigIndex j = start; j < end ; j ++ ) {
2100        CoinSimplexInt jColumn = indexColumnU[j];
2101        if (jColumn==pivotRow) {
2102          // swap
2103          numberInRowU[jRow]--;
2104          elementRowU[j]=elementRowU[end-1];
2105#if CONVERTROW
2106          convertRowToColumn[j]=convertRowToColumn[end-1];
2107#endif
2108          indexColumnU[j]=indexColumnU[end-1];
2109          break;
2110        }
2111      }
2112    }   
2113  }
2114#endif   
2115  //zero out pivot Row (before or after?)
2116  //add to R
2117  CoinBigIndex * COIN_RESTRICT startColumnR = startColumnRAddress_;
2118  CoinBigIndex putR = lengthR_;
2119  CoinSimplexInt number = numberR_;
2120 
2121  startColumnR[number] = putR;  //for luck and first time
2122  number++;
2123  //assert (startColumnR+number-firstCountAddress_<( CoinMax(5*numberRows_,2*numberRows_+2*maximumPivots_)+2));
2124  startColumnR[number] = putR + numberNonZero;
2125  numberR_ = number;
2126  lengthR_ = putR + numberNonZero;
2127  totalElements_ += numberNonZero;
2128
2129#if ABC_SMALL<2
2130  CoinSimplexInt *  COIN_RESTRICT nextRow=NULL;
2131  CoinSimplexInt *  COIN_RESTRICT lastRow=NULL;
2132  CoinSimplexInt next;
2133  CoinSimplexInt last;
2134  if (gotUCopy()) {
2135    //take out row
2136    nextRow = nextRowAddress_;
2137    lastRow = lastRowAddress_;
2138    next = nextRow[pivotRow];
2139    last = lastRow[pivotRow];
2140   
2141    nextRow[last] = next;
2142    lastRow[next] = last;
2143    numberInRowU[pivotRow]=0;
2144  }
2145#endif 
2146  //put in pivot
2147  //add row counts
2148 
2149  int n = 0;
2150  for (CoinSimplexInt i = 0; i < numberInColumnU2; i++ ) {
2151    CoinSimplexInt iRow = indexU2[i];
2152    CoinFactorizationDouble value=elementU2[i];
2153    assert (value);
2154    if ( !TEST_LESS_THAN_TOLERANCE(value ) ) {
2155      if ( iRow != pivotRow ) {
2156        //modify by pivot
2157        CoinFactorizationDouble value2 = value*pivotValue;
2158        indexU2[n]=iRow;
2159        elementU2[n++] = value2;
2160#if ABC_SMALL<2
2161        if (gotUCopy()) {
2162          CoinSimplexInt next = nextRow[iRow];
2163          CoinSimplexInt iNumberInRow = numberInRowU[iRow];
2164          CoinBigIndex space;
2165          CoinBigIndex put = startRowU[iRow] + iNumberInRow;
2166          space = startRowU[next] - put;
2167          if ( space <= 0 ) {
2168            getRowSpaceIterate ( iRow, iNumberInRow + 4 );
2169            put = startRowU[iRow] + iNumberInRow;
2170          } else if (next==numberRows_) {
2171            lastEntryByRowU_=put+1;
2172          }     
2173          indexColumnU[put] = pivotRow;
2174#if CONVERTROW
2175          convertRowToColumn[put] = i;
2176#endif
2177          elementRowU[put] = value2;
2178          numberInRowU[iRow] = iNumberInRow + 1;
2179        }
2180#endif
2181      } else {
2182        //swap and save
2183        indexU2[i]=indexU2[numberInColumnU2-1];
2184        elementU2[i] = elementU2[numberInColumnU2-1];
2185        numberInColumnU2--;
2186        i--;
2187      }
2188    } else {
2189      // swap
2190      indexU2[i]=indexU2[numberInColumnU2-1];
2191      elementU2[i] = elementU2[numberInColumnU2-1];
2192      numberInColumnU2--;
2193      i--;
2194    }       
2195  }       
2196  numberInColumnU2=n;
2197  //do permute
2198  permuteAddress_[numberRowsExtra_] = pivotRow;
2199  //and for safety
2200  permuteAddress_[numberRowsExtra_ + 1] = 0;
2201
2202  //pivotColumnAddress_[pivotRow] = numberRowsExtra_;
2203 
2204  numberU_++;
2205
2206#if ABC_SMALL<2
2207  if (gotUCopy()) {
2208    //in at end
2209    last = lastRow[numberRows_];
2210    nextRow[last] = pivotRow; //numberRowsExtra_;
2211    lastRow[numberRows_] = pivotRow; //numberRowsExtra_;
2212    lastRow[pivotRow] = last;
2213    nextRow[pivotRow] = numberRows_;
2214    startRowU[pivotRow] = lastEntryByRowU_;
2215  }
2216#endif
2217#ifndef ABC_USE_FUNCTION_POINTERS
2218  numberInColumn[pivotRow]=numberInColumnU2;
2219#else
2220#if ABC_USE_FUNCTION_POINTERS
2221  if (numberInColumnU2<9) {
2222    scatter[pivotRow].functionPointer=AbcScatterLowSubtract[numberInColumnU2];
2223  } else {
2224    scatter[pivotRow].functionPointer=AbcScatterHighSubtract[numberInColumnU2&3];
2225  }
2226#endif
2227  assert(scatter[pivotRow].offset==lastEntryByColumnUPlus_);
2228  //CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+lastEntryByColumnUPlus_;
2229  //CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+numberInColumnU2);
2230  if (numberInColumnU2<scatter[numberRows_].number) {
2231    int offset=2*(scatter[numberRows_].number-numberInColumnU2);
2232    indexU2 -= offset;
2233    //assert(indexU2-reinterpret_cast<int*>(area)>=0);
2234    for (int i=0;i<numberInColumnU2;i++)
2235      indexU2[i]=indexU2[i+offset];
2236  }
2237  scatter[pivotRow].number=numberInColumnU2;
2238  //CoinAbcMemcpy(indices,indexU2,numberInColumnU2);
2239  //CoinAbcMemcpy(area,elementU2,numberInColumnU2);
2240  lastEntryByColumnUPlus_ += (3*numberInColumnU2+1)>>1;
2241#endif
2242  totalElements_ += numberInColumnU2;
2243  lengthU_ += numberInColumnU2;
2244#if ABC_SMALL<2
2245  CoinSimplexInt * COIN_RESTRICT nextColumn = nextColumnAddress_;
2246  CoinSimplexInt * COIN_RESTRICT lastColumn = lastColumnAddress_;
2247  if (gotRCopy()) {
2248    //column in at beginning (as empty)
2249    next = nextColumn[maximumRowsExtra_];
2250    lastColumn[next] = numberRowsExtra_;
2251    nextColumn[maximumRowsExtra_] = numberRowsExtra_;
2252    nextColumn[numberRowsExtra_] = next;
2253    lastColumn[numberRowsExtra_] = maximumRowsExtra_;
2254  }
2255#endif
2256#if 0
2257  {
2258    int k=-1;
2259    for (int i=0;i<numberRows_;i++) {
2260      k=CoinMax(k,startRowU[i]+numberInRowU[i]);
2261    }
2262    assert (k==lastEntryByRowU_);
2263  }
2264#endif
2265  //pivotRegion[numberRowsExtra_] = pivotValue;
2266  pivotRegion[pivotRow] = pivotValue;
2267#ifndef ABC_USE_FUNCTION_POINTERS
2268  maximumU_ = CoinMax(maximumU_,startU+numberInColumnU2);
2269#else
2270  maximumU_ = CoinMax(maximumU_,lastEntryByColumnUPlus_);
2271#endif
2272  permuteLookup[pivotRow]=numberRowsExtra_;
2273  permuteLookup[numberRowsExtra_]=pivotRow;
2274  //pivotColumnAddress_[pivotRow]=numberRowsExtra_;
2275  //pivotColumnAddress_[numberRowsExtra_]=pivotRow;
2276  assert (pivotRow<numberRows_);
2277  // out of chain
2278  CoinSimplexInt iLast=pivotLinkedBackwardsAddress_[pivotRow];
2279  CoinSimplexInt iNext=pivotLinkedForwardsAddress_[pivotRow];
2280  assert (pivotRow==pivotLinkedForwardsAddress_[iLast]);
2281  assert (pivotRow==pivotLinkedBackwardsAddress_[iNext]);
2282  pivotLinkedForwardsAddress_[iLast]=iNext;
2283  pivotLinkedBackwardsAddress_[iNext]=iLast;
2284  if (pivotRow==lastSlack_) {
2285    lastSlack_ = iLast;
2286  }
2287  iLast=pivotLinkedBackwardsAddress_[numberRows_];
2288  assert (numberRows_==pivotLinkedForwardsAddress_[iLast]);
2289  pivotLinkedForwardsAddress_[iLast]=pivotRow;
2290  pivotLinkedBackwardsAddress_[pivotRow]=iLast;
2291  pivotLinkedBackwardsAddress_[numberRows_]=pivotRow;
2292  pivotLinkedForwardsAddress_[pivotRow]=numberRows_;
2293  assert (numberRows_>pivotLinkedForwardsAddress_[-1]);
2294  assert (pivotLinkedBackwardsAddress_[numberRows_]>=0);
2295  numberRowsExtra_++;
2296  numberGoodU_++;
2297  numberPivots_++;
2298  if ( numberRowsExtra_ > numberRows_ + 50 ) {
2299    CoinBigIndex extra = factorElements_ >> 1;
2300   
2301    if ( numberRowsExtra_ > numberRows_ + 100 + numberRows_ / 500 ) {
2302      if ( extra < 2 * numberRows_ ) {
2303        extra = 2 * numberRows_;
2304      }       
2305    } else {
2306      if ( extra < 5 * numberRows_ ) {
2307        extra = 5 * numberRows_;
2308      }       
2309    }       
2310    CoinBigIndex added = totalElements_ - factorElements_;
2311   
2312    if ( added > extra && added > ( factorElements_ ) << 1 && !status
2313         && 3*totalElements_ > 2*(lengthAreaU_+lengthAreaL_)) {
2314      status = 3;
2315      if ( messageLevel_ & 4 ) {
2316        std::cout << "Factorization has "<< totalElements_
2317                  << ", basis had " << factorElements_ <<std::endl;
2318      }
2319    }       
2320  }
2321#if ABC_SMALL<2
2322  if (gotRCopy()&&status<2) {
2323    //if (numberInColumnPlus&&status<2) {
2324    // we are going to put another copy of R in R
2325    CoinFactorizationDouble * COIN_RESTRICT elementRR = elementRAddress_ + lengthAreaR_;
2326    CoinSimplexInt * COIN_RESTRICT indexRowRR = indexRowRAddress_ + lengthAreaR_;
2327    CoinBigIndex * COIN_RESTRICT startRR = startColumnRAddress_+maximumPivots_+1;
2328    CoinSimplexInt pivotRow = numberRowsExtra_-1;
2329    for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
2330      CoinSimplexInt jRow = regionIndex[i];
2331      CoinFactorizationDouble value=region[jRow];
2332      elementR[putR] = value;
2333      indexRowR[putR] = jRow;
2334      putR++;
2335      region[jRow]=0.0;
2336      int iRow = permuteLookup[jRow];
2337      next = nextColumn[iRow];
2338      CoinBigIndex space;
2339      if (next!=maximumRowsExtra_)
2340        space = startRR[next]-startRR[iRow];
2341      else
2342        space = lengthAreaR_-startRR[iRow];
2343      CoinSimplexInt numberInR = numberInColumnPlus[iRow];
2344      if (space>numberInR) {
2345        // there is space
2346        CoinBigIndex  put=startRR[iRow]+numberInR;
2347        numberInColumnPlus[iRow]=numberInR+1;
2348        indexRowRR[put]=pivotRow;
2349        elementRR[put]=value;
2350        //add 4 for luck
2351        if (next==maximumRowsExtra_)
2352          startRR[maximumRowsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaR_);
2353      } else {
2354        // no space - do we shuffle?
2355        if (!getColumnSpaceIterateR(iRow,value,pivotRow)) {
2356          // printf("Need more space for R\n");
2357          numberInColumnPlus_.conditionalDelete();
2358          numberInColumnPlusAddress_=NULL;
2359          setNoGotRCopy();
2360          regionSparse->clear();
2361#if ABC_SMALL<0
2362          abort();
2363          status=3;
2364#endif
2365          break;
2366        }
2367      }
2368    }       
2369  } else {
2370#endif
2371    for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
2372      CoinSimplexInt iRow = regionIndex[i];
2373      elementR[putR] = region[iRow];
2374      region[iRow]=0.0;
2375      indexRowR[putR] = iRow;
2376      putR++;
2377    }       
2378#if ABC_SMALL<2
2379  }
2380#endif
2381#ifdef ABC_USE_FUNCTION_POINTERS
2382#if 0
2383  {
2384    CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
2385    CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
2386    CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
2387    CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
2388    CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
2389    scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
2390    extern scatterUpdate AbcScatterLowSubtract[9];
2391    extern scatterUpdate AbcScatterHighSubtract[4];
2392    printf("=============================================start\n");
2393    for (int iRow=0;iRow<numberRows_;iRow++) {
2394      int number=scatter[iRow].number;
2395      CoinBigIndex offset = scatter[iRow].offset;
2396      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+offset;
2397      CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+number);
2398      CoinSimplexInt * COIN_RESTRICT indices2 = indexRow+startColumnU[iRow];
2399      CoinFactorizationDouble * COIN_RESTRICT area2 = element+startColumnU[iRow];
2400      //assert (number==numberInColumn[iRow]);
2401#if ABC_USE_FUNCTION_POINTERS
2402      if (number<9)
2403        assert (scatter[iRow].functionPointer==AbcScatterLowSubtract[number]);
2404      else
2405        assert (scatter[iRow].functionPointer==AbcScatterHighSubtract[number&3]);
2406#endif
2407      if (number)
2408        printf("Row %d %d els\n",iRow,number);
2409      int n=0;
2410      for (int i=0;i<number;i++) {
2411        printf("(%d,%g) ",indices[i],area[i]);
2412        n++;
2413        if (n==10) {
2414          n=0;
2415          printf("\n");
2416        }
2417      }
2418      if (n)
2419        printf("\n");
2420#if 0
2421      for (int i=0;i<number;i++) {
2422        assert (indices[i]==indices2[i]);
2423        assert (area[i]==area2[i]);
2424      }
2425#endif
2426    }
2427    printf("=============================================end\n");
2428  }
2429#endif
2430#endif
2431  regionSparse->setNumElements(0);
2432  fromLongArray(regionSparse);
2433#if 0
2434  static int arrayF[6]={1,1,1,1,1,1};
2435  static int arrayB[6]={1,1,1,1,1,1};
2436  static int arrayF1[6]={1,1,1,1,1,1};
2437  static int arrayB1[6]={1,1,1,1,1,1};
2438  static int itry=0;
2439  {
2440    const CoinSimplexInt *  COIN_RESTRICT back = firstCountAddress_+numberRows_;
2441    const CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_; // no
2442    const CoinSimplexInt *  COIN_RESTRICT permute = permuteAddress_;
2443    const CoinSimplexInt *  COIN_RESTRICT pivotLOrder = this->pivotLOrderAddress_;
2444    const CoinSimplexInt *  COIN_RESTRICT pivotLinkedForwards = this->pivotLinkedForwardsAddress_;//no
2445    const CoinSimplexInt *  COIN_RESTRICT pivotLinkedBackwards = this->pivotLinkedBackwardsAddress_;//no
2446    itry++;
2447    int testit[6];
2448    for (int i=0;i<6;i++)
2449      testit[i]=1;
2450    for (int i=0;i<numberRows_;i++) {
2451      if (back[i]!=i)
2452        testit[0]=0;
2453      if (pivotColumn[i]!=i)
2454        testit[1]=0;
2455      if (permute[i]!=i)
2456        testit[2]=0;
2457      if (pivotLOrder[i]!=i)
2458        testit[3]=0;
2459      if (pivotLinkedForwards[i]!=i)
2460        testit[4]=0;
2461      if (pivotLinkedBackwards[i]!=i)
2462        testit[5]=0;
2463    }
2464    for (int i=0;i<6;i++) {
2465      if (!testit[i])
2466        arrayF[i]=0;
2467    }
2468    for (int i=0;i<6;i++)
2469      testit[i]=1;
2470    for (int i=0;i<numberRows_;i++) {
2471      if (back[i]!=(numberRows_-1-i))
2472        testit[0]=0;
2473      if (pivotColumn[i]!=(numberRows_-1-i))
2474        testit[1]=0;
2475      if (permute[i]!=(numberRows_-1-i))
2476        testit[2]=0;
2477      if (pivotLOrder[i]!=(numberRows_-1-i))
2478        testit[3]=0;
2479      if (pivotLinkedForwards[i]!=(numberRows_-1-i))
2480        testit[4]=0;
2481      if (pivotLinkedBackwards[i]!=(numberRows_-1-i))
2482        testit[5]=0;
2483    }
2484    for (int i=0;i<6;i++) {
2485      if (!testit[i])
2486        arrayB[i]=0;
2487    }
2488    for (int i=0;i<6;i++)
2489      testit[i]=1;
2490    for (int i=0;i<numberRows_-1;i++) {
2491      if (back[i]!=(i+1))
2492        testit[0]=0;
2493      if (pivotColumn[i]!=(i+1))
2494        testit[1]=0;
2495      if (permute[i]!=(i+1))
2496        testit[2]=0;
2497      if (pivotLOrder[i]!=(i+1))
2498        testit[3]=0;
2499      if (pivotLinkedForwards[i]!=(i+1))
2500        testit[4]=0;
2501      if (pivotLinkedBackwards[i]!=(i+1))
2502        testit[5]=0;
2503    }
2504    for (int i=0;i<6;i++) {
2505      if (!testit[i])
2506        arrayF1[i]=0;
2507    }
2508    for (int i=0;i<6;i++)
2509      testit[i]=1;
2510    for (int i=0;i<numberRows_-1;i++) {
2511      if (back[i]!=pivotLinkedBackwards[i-1])
2512        testit[0]=0;
2513      if (pivotColumn[i]!=(numberRows_-i-2))
2514        testit[1]=0;
2515      if (permute[i]!=(numberRows_-i-2))
2516        testit[2]=0;
2517      if (pivotLOrder[i]!=(numberRows_-i-2))
2518        testit[3]=0;
2519      if (pivotLinkedForwards[i]!=(numberRows_-i-2))
2520        testit[4]=0;
2521      if (pivotLinkedBackwards[i]!=i-1)
2522        testit[5]=0;
2523    }
2524    for (int i=0;i<6;i++) {
2525      if (!testit[i])
2526        arrayB1[i]=0;
2527    }
2528    if ((itry%1000)==0) {
2529      for (int i=0;i<6;i++) {
2530        if (arrayF[i])
2531          printf("arrayF[%d] still active\n",i);
2532        if (arrayB[i])
2533          printf("arrayB[%d] still active\n",i);
2534        if (arrayF1[i])
2535          printf("arrayF1[%d] still active\n",i);
2536        if (arrayB1[i])
2537          printf("arrayB1[%d] still active\n",i);
2538      }
2539    }
2540  }
2541#endif
2542}
2543/* Replaces one Column to basis,
2544   partial update in vector */
2545void 
2546CoinAbcTypeFactorization::replaceColumnPart3 ( const AbcSimplex * /*model*/,
2547                                          CoinIndexedVector * regionSparse,
2548                                               CoinIndexedVector * /*tableauColumn*/,
2549                                               CoinIndexedVector * partialUpdate,
2550                                          int pivotRow,
2551#ifdef ABC_LONG_FACTORIZATION
2552                                               long
2553#endif
2554                                          double alpha )
2555{
2556  assert (numberU_<=numberRowsExtra_);
2557#ifndef ABC_USE_FUNCTION_POINTERS
2558  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
2559  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
2560  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
2561#endif
2562  CoinSimplexInt * COIN_RESTRICT indexRowR = indexRowRAddress_;
2563#if ABC_SMALL<2
2564  CoinSimplexInt * COIN_RESTRICT numberInRowU = numberInRowAddress_;
2565#endif
2566#if ABC_SMALL<2
2567  CoinSimplexInt * COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
2568#endif
2569  CoinSimplexInt *  COIN_RESTRICT permuteLookup = pivotColumnAddress_;
2570  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
2571 
2572  CoinSimplexInt numberInColumnU2 = partialUpdate->getNumElements();
2573  //take out old pivot column
2574 
2575#ifndef ABC_USE_FUNCTION_POINTERS
2576  totalElements_ -= numberInColumn[pivotRow];
2577  CoinBigIndex startU = startColumnU[numberRows_];
2578  startColumnU[pivotRow]=startU;
2579  CoinSimplexInt * COIN_RESTRICT indexU2 = &indexRowUAddress_[startU];
2580  CoinFactorizationDouble * COIN_RESTRICT elementU2 = &elementUAddress_[startU];
2581  numberInColumn[numberRows_]=numberInColumnU2;
2582#else
2583  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
2584#if ABC_USE_FUNCTION_POINTERS
2585  extern scatterUpdate AbcScatterLowSubtract[9];
2586  extern scatterUpdate AbcScatterHighSubtract[4];
2587#endif
2588  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
2589  CoinBigIndex startU = lastEntryByColumnUPlus_;
2590  totalElements_ -= scatter[pivotRow].number;
2591  CoinBigIndex saveStart = scatter[pivotRow].offset;
2592  CoinBigIndex saveEnd = scatter[pivotRow].number;
2593  scatter[pivotRow].offset=startU;
2594
2595  scatter[numberRows_].number=numberInColumnU2;
2596  CoinFactorizationDouble * COIN_RESTRICT elementU2 = elementUColumnPlus+startU;
2597  CoinSimplexInt * COIN_RESTRICT indexU2 = reinterpret_cast<CoinSimplexInt *>(elementU2+numberInColumnU2);
2598#endif
2599  // move
2600  CoinSimplexInt * COIN_RESTRICT indexU2Save = partialUpdate->getIndices();
2601  CoinFactorizationDouble * COIN_RESTRICT elementU2Save = denseVector(partialUpdate);
2602  memcpy(indexU2,indexU2Save,numberInColumnU2*sizeof(CoinSimplexInt));
2603  for (int i=0;i<numberInColumnU2;i++) {
2604    int iRow=indexU2[i];
2605    elementU2[i]=elementU2Save[iRow];
2606    elementU2Save[iRow]=0.0;
2607  }
2608  //memcpy(elementU2,elementU2Save,numberInColumnU2*sizeof(double));
2609  //memset(elementU2Save,0,numberInColumnU2*sizeof(double));
2610  partialUpdate->setNumElements(0);
2611  //get entries in row (pivot not stored)
2612  CoinSimplexInt numberNonZero = 0;
2613#if ABC_SMALL<2
2614  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
2615#endif
2616#if CONVERTROW
2617  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
2618#if CONVERTROW>2
2619  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
2620#endif
2621#endif
2622#if ABC_SMALL<2
2623  CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
2624  CoinBigIndex * COIN_RESTRICT startRowU = startRowUAddress_;
2625#endif
2626  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
2627  CoinFactorizationDouble pivotValue = 1.0;
2628  CoinSimplexInt status=0;
2629#ifndef ABC_USE_FUNCTION_POINTERS
2630  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
2631#endif
2632  CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_;
2633
2634#ifndef NDEBUG
2635#if CONVERTROW
2636#define CONVERTDEBUG
2637#endif
2638#endif
2639  int nInRow=9999999; // ? what if ABC_SMALL==0 or ABC_SMALL==1
2640#if ABC_SMALL<2
2641  CoinBigIndex start=0;
2642  CoinBigIndex end=0;
2643  if (gotUCopy()) {
2644    start=startRowU[pivotRow];
2645    nInRow=numberInRowU[pivotRow];
2646    end= start + nInRow;
2647  }
2648#endif
2649  CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
2650  numberNonZero = regionSparse->getNumElements (  );
2651  //CoinFactorizationDouble saveFromU = 0.0;
2652
2653  instrument_do("CoinAbcFactorizationReplaceColumn",2*numberRows_+4*numberInColumnU2);
2654  pivotValue = 1.0 / alpha;
2655#if ABC_SMALL<2
2656  if (gotUCopy()) {
2657    for (CoinBigIndex i = start; i < end ; i ++ ) {
2658      CoinSimplexInt jColumn = indexColumnU[i];
2659#ifndef ABC_USE_FUNCTION_POINTERS
2660      CoinBigIndex startColumn = startColumnU[jColumn];
2661#if CONVERTROW
2662      CoinBigIndex j = convertRowToColumn[i];
2663#else
2664      CoinBigIndex j;
2665      int number = numberInColumn[jColumn];
2666      for (j=0;j<number;j++) {
2667        if (indexRowU[j+startColumn]==pivotRow)
2668          break;
2669      }
2670      assert (j<number);
2671      //assert (j==convertRowToColumn[i]); // temp
2672#endif
2673#ifdef CONVERTDEBUG
2674      assert (fabs(elementU[j+startColumn]-elementRowU[i])<1.0e-4);
2675#endif     
2676#else
2677      int number = scatter[jColumn].number;
2678      CoinSimplexInt k = number-1;
2679      scatter[jColumn].number=k;
2680      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter[jColumn].offset;
2681      CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+k+1);
2682#if CONVERTROW
2683      CoinSimplexInt j = convertRowToColumn[i];
2684#else
2685      CoinSimplexInt j;
2686      for (j=0;j<number;j++) {
2687        if (indices[j]==pivotRow)
2688          break;
2689      }
2690      assert (j<number);
2691      //assert (j==convertRowToColumn[i]); // temp
2692#endif
2693#ifdef CONVERTDEBUG
2694      assert (fabs(area[j]-elementRowU[i])<1.0e-4);
2695#endif     
2696#endif
2697      // swap
2698#ifndef ABC_USE_FUNCTION_POINTERS
2699      CoinSimplexInt k = numberInColumn[jColumn]-1;
2700      numberInColumn[jColumn]=k;
2701      CoinBigIndex k2 = k+startColumn;
2702      int kRow2=indexRowU[k2];
2703      indexRowU[j+startColumn]=kRow2;
2704      CoinFactorizationDouble value2=elementU[k2];
2705      elementU[j+startColumn]=value2;
2706#else
2707      int kRow2=indices[k];
2708      CoinFactorizationDouble value2=area[k];
2709#if ABC_USE_FUNCTION_POINTERS
2710      if (k<9) {
2711        scatter[jColumn].functionPointer=AbcScatterLowSubtract[k];
2712      } else {
2713        scatter[jColumn].functionPointer=AbcScatterHighSubtract[k&3];
2714      }
2715#endif
2716      // later more complicated swap to avoid move (don't think we can!)
2717      indices[j]=kRow2;
2718      area[j]=value2;
2719      // move
2720      indices-=2;
2721      for (int i=0;i<k;i++)
2722        indices[i]=indices[i+2];
2723#endif
2724      // move in row copy (slow - as other remark (but shouldn't need to))
2725      CoinBigIndex start2 = startRowU[kRow2];
2726      int n=numberInRowU[kRow2];
2727      CoinBigIndex end2 = start2 + n;
2728#ifndef NDEBUG
2729      bool found=false;
2730#endif
2731      for (CoinBigIndex i2 = start2; i2 < end2 ; i2 ++ ) {
2732        CoinSimplexInt iColumn2 = indexColumnU[i2];
2733        if (jColumn==iColumn2) {
2734#if CONVERTROW
2735          convertRowToColumn[i2] = j;
2736#endif
2737#ifndef NDEBUG
2738          found=true;
2739#endif
2740          break;
2741        }
2742      }
2743#ifndef NDEBUG
2744      assert(found);
2745#endif
2746    }
2747  } else {
2748#endif
2749#if ABC_SMALL>=0
2750    assert(nInRow!=9999999); // ? what if ABC_SMALL==0 or ABC_SMALL==1
2751    // no row copy
2752    CoinBigIndex COIN_RESTRICT * deletedPosition = reinterpret_cast<CoinBigIndex *>(elementR+lengthR_);
2753    CoinSimplexInt COIN_RESTRICT * deletedColumns = reinterpret_cast<CoinSimplexInt *>(indexRowR+lengthR_);
2754    for (CoinSimplexInt i = 0; i < nInRow ; i ++ ) {
2755      CoinBigIndex j = deletedPosition[i];
2756      CoinSimplexInt jColumn = deletedColumns[i];
2757      // swap
2758      CoinSimplexInt k = numberInColumn[jColumn]-1;
2759      numberInColumn[jColumn]=k;
2760      CoinBigIndex k2 = k+startColumnU[jColumn];
2761      int kRow2=indexRowU[k2];
2762      indexRowU[j]=kRow2;
2763      CoinFactorizationDouble value2=elementU[k2];
2764      elementU[j]=value2;
2765    }
2766#endif
2767#if ABC_SMALL<2
2768  }
2769#endif
2770 #if ABC_SMALL<2
2771  if (gotUCopy()) {
2772    // Now zero out column of U
2773    //take out old pivot column
2774#ifdef ABC_USE_FUNCTION_POINTERS
2775    CoinFactorizationDouble * COIN_RESTRICT elementU = elementUColumnPlus+saveStart;
2776    saveStart=0;
2777    CoinSimplexInt * COIN_RESTRICT indexRowU = reinterpret_cast<CoinSimplexInt *>(elementU+saveEnd);
2778#endif
2779    for (CoinBigIndex i = saveStart; i < saveEnd ; i ++ ) {
2780      //elementU[i] = 0.0;
2781      // If too slow then reverse meaning of convertRowToColumn and use
2782      int jRow=indexRowU[i];
2783      CoinBigIndex start = startRowU[jRow];
2784      CoinBigIndex end = start + numberInRowU[jRow];
2785      for (CoinBigIndex j = start; j < end ; j ++ ) {
2786        CoinSimplexInt jColumn = indexColumnU[j];
2787        if (jColumn==pivotRow) {
2788          // swap
2789          numberInRowU[jRow]--;
2790          elementRowU[j]=elementRowU[end-1];
2791#if CONVERTROW
2792          convertRowToColumn[j]=convertRowToColumn[end-1];
2793#endif
2794          indexColumnU[j]=indexColumnU[end-1];
2795          break;
2796        }
2797      }
2798    }   
2799  }
2800#endif   
2801  //zero out pivot Row (before or after?)
2802  //add to R
2803  CoinBigIndex * COIN_RESTRICT startColumnR = startColumnRAddress_;
2804  CoinBigIndex putR = lengthR_;
2805  CoinSimplexInt number = numberR_;
2806 
2807  startColumnR[number] = putR;  //for luck and first time
2808  number++;
2809  //assert (startColumnR+number-firstCountAddress_<( CoinMax(5*numberRows_,2*numberRows_+2*maximumPivots_)+2));
2810  startColumnR[number] = putR + numberNonZero;
2811  numberR_ = number;
2812  lengthR_ = putR + numberNonZero;
2813  totalElements_ += numberNonZero;
2814
2815#if ABC_SMALL<2
2816  CoinSimplexInt *  COIN_RESTRICT nextRow=NULL;
2817  CoinSimplexInt *  COIN_RESTRICT lastRow=NULL;
2818  CoinSimplexInt next;
2819  CoinSimplexInt last;
2820  if (gotUCopy()) {
2821    //take out row
2822    nextRow = nextRowAddress_;
2823    lastRow = lastRowAddress_;
2824    next = nextRow[pivotRow];
2825    last = lastRow[pivotRow];
2826   
2827    nextRow[last] = next;
2828    lastRow[next] = last;
2829    numberInRowU[pivotRow]=0;
2830  }
2831#endif 
2832  //put in pivot
2833  //add row counts
2834 
2835  int n = 0;
2836  for (CoinSimplexInt i = 0; i < numberInColumnU2; i++ ) {
2837    CoinSimplexInt iRow = indexU2[i];
2838    CoinFactorizationDouble value=elementU2[i];
2839    assert (value);
2840    if ( !TEST_LESS_THAN_TOLERANCE(value ) ) {
2841      if ( iRow != pivotRow ) {
2842        //modify by pivot
2843        CoinFactorizationDouble value2 = value*pivotValue;
2844        indexU2[n]=iRow;
2845        elementU2[n++] = value2;
2846#if ABC_SMALL<2
2847        if (gotUCopy()) {
2848          CoinSimplexInt next = nextRow[iRow];
2849          CoinSimplexInt iNumberInRow = numberInRowU[iRow];
2850          CoinBigIndex space;
2851          CoinBigIndex put = startRowU[iRow] + iNumberInRow;
2852          space = startRowU[next] - put;
2853          if ( space <= 0 ) {
2854            getRowSpaceIterate ( iRow, iNumberInRow + 4 );
2855            put = startRowU[iRow] + iNumberInRow;
2856          } else if (next==numberRows_) {
2857            lastEntryByRowU_=put+1;
2858          }     
2859          indexColumnU[put] = pivotRow;
2860#if CONVERTROW
2861          convertRowToColumn[put] = i;
2862#endif
2863          elementRowU[put] = value2;
2864          numberInRowU[iRow] = iNumberInRow + 1;
2865        }
2866#endif
2867      } else {
2868        //swap and save
2869        indexU2[i]=indexU2[numberInColumnU2-1];
2870        elementU2[i] = elementU2[numberInColumnU2-1];
2871        numberInColumnU2--;
2872        i--;
2873      }
2874    } else {
2875      // swap
2876      indexU2[i]=indexU2[numberInColumnU2-1];
2877      elementU2[i] = elementU2[numberInColumnU2-1];
2878      numberInColumnU2--;
2879      i--;
2880    }       
2881  }       
2882  numberInColumnU2=n;
2883  //do permute
2884  permuteAddress_[numberRowsExtra_] = pivotRow;
2885  //and for safety
2886  permuteAddress_[numberRowsExtra_ + 1] = 0;
2887
2888  //pivotColumnAddress_[pivotRow] = numberRowsExtra_;
2889 
2890  numberU_++;
2891
2892#if ABC_SMALL<2
2893  if (gotUCopy()) {
2894    //in at end
2895    last = lastRow[numberRows_];
2896    nextRow[last] = pivotRow; //numberRowsExtra_;
2897    lastRow[numberRows_] = pivotRow; //numberRowsExtra_;
2898    lastRow[pivotRow] = last;
2899    nextRow[pivotRow] = numberRows_;
2900    startRowU[pivotRow] = lastEntryByRowU_;
2901  }
2902#endif
2903#ifndef ABC_USE_FUNCTION_POINTERS
2904  numberInColumn[pivotRow]=numberInColumnU2;
2905#else
2906#if ABC_USE_FUNCTION_POINTERS
2907  if (numberInColumnU2<9) {
2908    scatter[pivotRow].functionPointer=AbcScatterLowSubtract[numberInColumnU2];
2909  } else {
2910    scatter[pivotRow].functionPointer=AbcScatterHighSubtract[numberInColumnU2&3];
2911  }
2912#endif
2913  assert(scatter[pivotRow].offset==lastEntryByColumnUPlus_);
2914  //CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+lastEntryByColumnUPlus_;
2915  //CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+numberInColumnU2);
2916  if (numberInColumnU2<scatter[numberRows_].number) {
2917    int offset=2*(scatter[numberRows_].number-numberInColumnU2);
2918    indexU2 -= offset;
2919    //assert(indexU2-reinterpret_cast<int*>(area)>=0);
2920    for (int i=0;i<numberInColumnU2;i++)
2921      indexU2[i]=indexU2[i+offset];
2922  }
2923  scatter[pivotRow].number=numberInColumnU2;
2924  //CoinAbcMemcpy(indices,indexU2,numberInColumnU2);
2925  //CoinAbcMemcpy(area,elementU2,numberInColumnU2);
2926  lastEntryByColumnUPlus_ += (3*numberInColumnU2+1)>>1;
2927#endif
2928  totalElements_ += numberInColumnU2;
2929  lengthU_ += numberInColumnU2;
2930#if ABC_SMALL<2
2931  CoinSimplexInt * COIN_RESTRICT nextColumn = nextColumnAddress_;
2932  CoinSimplexInt * COIN_RESTRICT lastColumn = lastColumnAddress_;
2933  if (gotRCopy()) {
2934    //column in at beginning (as empty)
2935    next = nextColumn[maximumRowsExtra_];
2936    lastColumn[next] = numberRowsExtra_;
2937    nextColumn[maximumRowsExtra_] = numberRowsExtra_;
2938    nextColumn[numberRowsExtra_] = next;
2939    lastColumn[numberRowsExtra_] = maximumRowsExtra_;
2940  }
2941#endif
2942#if 0
2943  {
2944    int k=-1;
2945    for (int i=0;i<numberRows_;i++) {
2946      k=CoinMax(k,startRowU[i]+numberInRowU[i]);
2947    }
2948    assert (k==lastEntryByRowU_);
2949  }
2950#endif
2951  //pivotRegion[numberRowsExtra_] = pivotValue;
2952  pivotRegion[pivotRow] = pivotValue;
2953#ifndef ABC_USE_FUNCTION_POINTERS
2954  maximumU_ = CoinMax(maximumU_,startU+numberInColumnU2);
2955#else
2956  maximumU_ = CoinMax(maximumU_,lastEntryByColumnUPlus_);
2957#endif
2958  permuteLookup[pivotRow]=numberRowsExtra_;
2959  permuteLookup[numberRowsExtra_]=pivotRow;
2960  //pivotColumnAddress_[pivotRow]=numberRowsExtra_;
2961  //pivotColumnAddress_[numberRowsExtra_]=pivotRow;
2962  assert (pivotRow<numberRows_);
2963  // out of chain
2964  CoinSimplexInt iLast=pivotLinkedBackwardsAddress_[pivotRow];
2965  CoinSimplexInt iNext=pivotLinkedForwardsAddress_[pivotRow];
2966  assert (pivotRow==pivotLinkedForwardsAddress_[iLast]);
2967  assert (pivotRow==pivotLinkedBackwardsAddress_[iNext]);
2968  pivotLinkedForwardsAddress_[iLast]=iNext;
2969  pivotLinkedBackwardsAddress_[iNext]=iLast;
2970  if (pivotRow==lastSlack_) {
2971    lastSlack_ = iLast;
2972  }
2973  iLast=pivotLinkedBackwardsAddress_[numberRows_];
2974  assert (numberRows_==pivotLinkedForwardsAddress_[iLast]);
2975  pivotLinkedForwardsAddress_[iLast]=pivotRow;
2976  pivotLinkedBackwardsAddress_[pivotRow]=iLast;
2977  pivotLinkedBackwardsAddress_[numberRows_]=pivotRow;
2978  pivotLinkedForwardsAddress_[pivotRow]=numberRows_;
2979  assert (numberRows_>pivotLinkedForwardsAddress_[-1]);
2980  assert (pivotLinkedBackwardsAddress_[numberRows_]>=0);
2981  numberRowsExtra_++;
2982  numberGoodU_++;
2983  numberPivots_++;
2984  if ( numberRowsExtra_ > numberRows_ + 50 ) {
2985    CoinBigIndex extra = factorElements_ >> 1;
2986   
2987    if ( numberRowsExtra_ > numberRows_ + 100 + numberRows_ / 500 ) {
2988      if ( extra < 2 * numberRows_ ) {
2989        extra = 2 * numberRows_;
2990      }       
2991    } else {
2992      if ( extra < 5 * numberRows_ ) {
2993        extra = 5 * numberRows_;
2994      }       
2995    }       
2996    CoinBigIndex added = totalElements_ - factorElements_;
2997   
2998    if ( added > extra && added > ( factorElements_ ) << 1 && !status
2999         && 3*totalElements_ > 2*(lengthAreaU_+lengthAreaL_)) {
3000      status = 3;
3001      if ( messageLevel_ & 4 ) {
3002        std::cout << "Factorization has "<< totalElements_
3003                  << ", basis had " << factorElements_ <<std::endl;
3004      }
3005    }       
3006  }
3007#if ABC_SMALL<2
3008  if (gotRCopy()&&status<2) {
3009    //if (numberInColumnPlus&&status<2) {
3010    // we are going to put another copy of R in R
3011    CoinFactorizationDouble * COIN_RESTRICT elementRR = elementRAddress_ + lengthAreaR_;
3012    CoinSimplexInt * COIN_RESTRICT indexRowRR = indexRowRAddress_ + lengthAreaR_;
3013    CoinBigIndex * COIN_RESTRICT startRR = startColumnRAddress_+maximumPivots_+1;
3014    CoinSimplexInt pivotRow = numberRowsExtra_-1;
3015    for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
3016      CoinSimplexInt jRow = regionIndex[i];
3017      CoinFactorizationDouble value=region[jRow];
3018      elementR[putR] = value;
3019      indexRowR[putR] = jRow;
3020      putR++;
3021      region[jRow]=0.0;
3022      int iRow = permuteLookup[jRow];
3023      next = nextColumn[iRow];
3024      CoinBigIndex space;
3025      if (next!=maximumRowsExtra_)
3026        space = startRR[next]-startRR[iRow];
3027      else
3028        space = lengthAreaR_-startRR[iRow];
3029      CoinSimplexInt numberInR = numberInColumnPlus[iRow];
3030      if (space>numberInR) {
3031        // there is space
3032        CoinBigIndex  put=startRR[iRow]+numberInR;
3033        numberInColumnPlus[iRow]=numberInR+1;
3034        indexRowRR[put]=pivotRow;
3035        elementRR[put]=value;
3036        //add 4 for luck
3037        if (next==maximumRowsExtra_)
3038          startRR[maximumRowsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaR_);
3039      } else {
3040        // no space - do we shuffle?
3041        if (!getColumnSpaceIterateR(iRow,value,pivotRow)) {
3042          // printf("Need more space for R\n");
3043          numberInColumnPlus_.conditionalDelete();
3044          numberInColumnPlusAddress_=NULL;
3045          setNoGotRCopy();
3046          regionSparse->clear();
3047#if ABC_SMALL<0
3048          abort();
3049          status=3;
3050#endif
3051          break;
3052        }
3053      }
3054    }       
3055  } else {
3056#endif
3057    for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
3058      CoinSimplexInt iRow = regionIndex[i];
3059      elementR[putR] = region[iRow];
3060      region[iRow]=0.0;
3061      indexRowR[putR] = iRow;
3062      putR++;
3063    }       
3064#if ABC_SMALL<2
3065  }
3066#endif
3067#ifdef ABC_USE_FUNCTION_POINTERS
3068#if 0
3069  {
3070    CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
3071    CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
3072    CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
3073    CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
3074    CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
3075    scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
3076    extern scatterUpdate AbcScatterLowSubtract[9];
3077    extern scatterUpdate AbcScatterHighSubtract[4];
3078    printf("=============================================start\n");
3079    for (int iRow=0;iRow<numberRows_;iRow++) {
3080      int number=scatter[iRow].number;
3081      CoinBigIndex offset = scatter[iRow].offset;
3082      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+offset;
3083      CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+number);
3084      CoinSimplexInt * COIN_RESTRICT indices2 = indexRow+startColumnU[iRow];
3085      CoinFactorizationDouble * COIN_RESTRICT area2 = element+startColumnU[iRow];
3086      //assert (number==numberInColumn[iRow]);
3087#if ABC_USE_FUNCTION_POINTERS
3088      if (number<9)
3089        assert (scatter[iRow].functionPointer==AbcScatterLowSubtract[number]);
3090      else
3091        assert (scatter[iRow].functionPointer==AbcScatterHighSubtract[number&3]);
3092#endif
3093      if (number)
3094        printf("Row %d %d els\n",iRow,number);
3095      int n=0;
3096      for (int i=0;i<number;i++) {
3097        printf("(%d,%g) ",indices[i],area[i]);
3098        n++;
3099        if (n==10) {
3100          n=0;
3101          printf("\n");
3102        }
3103      }
3104      if (n)
3105        printf("\n");
3106#if 0
3107      for (int i=0;i<number;i++) {
3108        assert (indices[i]==indices2[i]);
3109        assert (area[i]==area2[i]);
3110      }
3111#endif
3112    }
3113    printf("=============================================end\n");
3114  }
3115#endif
3116#endif
3117  regionSparse->setNumElements(0);
3118}
3119#if ABC_SMALL>=0
3120/* Combines BtranU and store which elements are to be deleted
3121 */
3122int 
3123CoinAbcTypeFactorization::replaceColumnU ( CoinIndexedVector * regionSparse,
3124                                   CoinBigIndex * deletedPosition,
3125                                   CoinSimplexInt * deletedColumns,
3126                                   CoinSimplexInt pivotRow)
3127{
3128  instrument_start("CoinAbcFactorizationReplaceColumnU",numberRows_);
3129  CoinSimplexDouble * COIN_RESTRICT region = regionSparse->denseVector();
3130  int * COIN_RESTRICT regionIndex = regionSparse->getIndices();
3131  const CoinBigIndex COIN_RESTRICT *startColumn = startColumnUAddress_;
3132  const CoinSimplexInt COIN_RESTRICT *indexRow = indexRowUAddress_;
3133  CoinFactorizationDouble COIN_RESTRICT *element = elementUAddress_;
3134  int numberNonZero = 0;
3135  const CoinSimplexInt COIN_RESTRICT *numberInColumn = numberInColumnAddress_;
3136  int nPut=0;
3137  const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
3138  CoinSimplexInt jRow=pivotLinked[-1];
3139  while (jRow!=numberRows_) {
3140    assert (!region[jRow]);
3141    CoinFactorizationDouble pivotValue = 0.0;
3142    instrument_add(numberInColumn[jRow]);
3143    for (CoinBigIndex  j= startColumn[jRow] ; 
3144         j < startColumn[jRow]+numberInColumn[jRow]; j++ ) {
3145      int iRow = indexRow[j];
3146      CoinFactorizationDouble value = element[j];
3147      if (iRow!=pivotRow) {
3148        pivotValue -= value * region[iRow];
3149      } else {
3150        assert (!region[iRow]);
3151        pivotValue += value;
3152        deletedColumns[nPut]=jRow;
3153        deletedPosition[nPut++]=j;
3154      }
3155    }       
3156    if ( !TEST_LESS_THAN_TOLERANCE_REGISTER( pivotValue ) ) {
3157      regionIndex[numberNonZero++] = jRow;
3158      region[jRow] = pivotValue;
3159    } else {
3160      region[jRow] = 0;
3161    }
3162    jRow=pivotLinked[jRow];
3163  }
3164  regionSparse->setNumElements(numberNonZero);
3165  instrument_end();
3166  return nPut;
3167}
3168/* Updates part of column transpose (BTRANU) by column
3169   assumes index is sorted i.e. region is correct */
3170void 
3171CoinAbcTypeFactorization::updateColumnTransposeUByColumn ( CoinIndexedVector * regionSparse,
3172                                                   CoinSimplexInt smallestIndex) const
3173{
3174  instrument_start("CoinAbcFactorizationUpdateTransposeUByColumn",numberRows_);
3175  CoinSimplexDouble * COIN_RESTRICT region = regionSparse->denseVector();
3176  const CoinBigIndex COIN_RESTRICT *startColumn = startColumnUAddress_;
3177  const CoinSimplexInt COIN_RESTRICT *indexRow = indexRowUAddress_;
3178  const CoinFactorizationDouble COIN_RESTRICT *element = elementUAddress_;
3179#if ABC_SMALL<3
3180  int numberNonZero = 0;
3181  int * COIN_RESTRICT regionIndex = regionSparse->getIndices();
3182#endif
3183  const CoinSimplexInt COIN_RESTRICT *numberInColumn = numberInColumnAddress_;
3184  //const CoinFactorizationDouble COIN_RESTRICT *pivotRegion = pivotRegionAddress_;
3185  const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
3186  CoinSimplexInt jRow=smallestIndex;
3187  while (jRow!=numberRows_) {
3188    CoinFactorizationDouble pivotValue = region[jRow];
3189    CoinBigIndex start=startColumn[jRow];
3190    instrument_add(numberInColumn[jRow]);
3191#ifndef INLINE_IT2
3192    for (CoinBigIndex  j= start ; j < start+numberInColumn[jRow]; j++ ) {
3193      int iRow = indexRow[j];
3194      CoinFactorizationDouble value = element[j];
3195      pivotValue -= value * region[iRow];
3196    }
3197#else
3198    pivotValue+=CoinAbcGatherUpdate(numberInColumn[jRow],element+start,indexRow+start,region);
3199#endif       
3200#if ABC_SMALL<3
3201    if ( !TEST_LESS_THAN_TOLERANCE_REGISTER( pivotValue )  ) {
3202      regionIndex[numberNonZero++] = jRow;
3203      region[jRow] = pivotValue;
3204    } else {
3205      region[jRow] = 0;
3206    }
3207#else
3208      region[jRow] = pivotValue;
3209#endif       
3210    jRow=pivotLinked[jRow];
3211  }
3212#if ABC_SMALL<3
3213  regionSparse->setNumElements(numberNonZero);
3214#endif
3215  instrument_end();
3216}
3217#endif
3218//  updateColumnTranspose.  Updates one column transpose (BTRAN)
3219CoinSimplexInt
3220CoinAbcTypeFactorization::updateColumnTranspose ( CoinIndexedVector & regionSparse) const
3221{
3222  toLongArray(&regionSparse,0);
3223  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3224  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse.getIndices();
3225  CoinSimplexInt numberNonZero = regionSparse.getNumElements();
3226  const CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
3227
3228#ifndef ABC_ORDERED_FACTORIZATION
3229  // can I move this
3230#ifndef INLINE_IT3
3231  for (CoinSimplexInt i = 0; i < numberNonZero; i ++ ) {
3232    CoinSimplexInt iRow = regionIndex[i];
3233    CoinSimplexDouble value = region[iRow];
3234    value *= pivotRegion[iRow];
3235    region[iRow] = value;
3236  }
3237#else
3238  multiplyIndexed(numberNonZero,pivotRegion,
3239                  regionIndex,region);
3240#endif
3241#else
3242  // Permute in for Btran
3243  permuteInForBtranAndMultiply(regionSparse);
3244#endif
3245  //  ******* U
3246  // Apply pivot region - could be combined for speed
3247  // Can only combine/move inside vector for sparse
3248  CoinSimplexInt smallestIndex=pivotLinkedForwardsAddress_[-1];
3249#if ABC_SMALL<2
3250  // copy of code inside transposeU
3251  bool goSparse=false;
3252#else
3253#define goSparse false
3254#endif
3255#if ABC_SMALL<2
3256  // Guess at number at end
3257  if (gotUCopy()) {
3258    assert (btranAverageAfterU_);
3259    CoinSimplexInt newNumber = 
3260      static_cast<CoinSimplexInt> (numberNonZero*btranAverageAfterU_*twiddleBtranFactor1());
3261    if (newNumber< sparseThreshold_)
3262      goSparse = true;
3263  }
3264#endif
3265#if 0
3266  static int xxxxxx=0;
3267  xxxxxx++;
3268  if (xxxxxx%20==0) {
3269    int smallestK=-1;
3270    int largestK=-1;
3271    //const CoinBigIndex COIN_RESTRICT *startColumn = startColumnUAddress_;
3272    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
3273    const CoinSimplexInt COIN_RESTRICT *numberInColumn = numberInColumnAddress_;
3274    for (CoinSimplexInt k = numberRows_-1 ; k >= 0; k-- ) {
3275      CoinSimplexInt i=pivotRowForward[k];
3276      if (numberInColumn[i]) {
3277        if (largestK<0)
3278          largestK=k;
3279        smallestK=k;
3280      }
3281    }
3282    printf("UByCol %d slacks largestK %d smallestK %d\n",numberSlacks_,largestK,smallestK);
3283    //const CoinBigIndex * COIN_RESTRICT startRow = startRowUAddress_;
3284    const CoinSimplexInt * COIN_RESTRICT numberInRow = numberInRowAddress_;
3285    smallestK=-1;
3286    largestK=-1;
3287    for (CoinSimplexInt k = numberRows_-1 ; k >= 0; k-- ) {
3288      CoinSimplexInt i=pivotRowForward[k];
3289      if (numberInRow[i]) {
3290        if (largestK<0)
3291          largestK=k;
3292        smallestK=k;
3293      }
3294    }
3295    printf("and ByRow largestK %d smallestK %d\n",largestK,smallestK);
3296  }
3297#endif
3298  if (numberNonZero<40&&(numberNonZero<<4)<numberRows_&&!goSparse) {
3299    CoinSimplexInt *  COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
3300    CoinSimplexInt smallest=numberRowsExtra_;
3301    for (CoinSimplexInt j = 0; j < numberNonZero; j++ ) {
3302      CoinSimplexInt iRow = regionIndex[j];
3303      if (pivotRowForward[iRow]<smallest) {
3304        smallest=pivotRowForward[iRow];
3305        smallestIndex=iRow;
3306      }
3307    }
3308  }
3309  updateColumnTransposeU ( &regionSparse,smallestIndex
3310#if ABC_SMALL<2
3311                  ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
3312#endif
3313#if ABC_PARALLEL
3314                    ,0
3315#endif
3316                           );
3317  //row bits here
3318  updateColumnTransposeR ( &regionSparse
3319#if ABC_SMALL<2
3320                  ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
3321#endif
3322                           );
3323  //  ******* L
3324  updateColumnTransposeL ( &regionSparse
3325#if ABC_SMALL<2
3326                  ,reinterpret_cast<CoinAbcStatistics &>(btranCountInput_)
3327#endif
3328#if ABC_PARALLEL
3329                                      ,0
3330#endif
3331                           );
3332#if ABC_SMALL<3
3333#ifdef ABC_ORDERED_FACTORIZATION
3334  // Permute out for Btran
3335  permuteOutForBtran(regionSparse);
3336#endif
3337#if ABC_DEBUG
3338  regionSparse.checkClean();
3339#endif
3340  numberNonZero = regionSparse.getNumElements (  );
3341#else
3342  numberNonZero=0;
3343  for (CoinSimplexInt i=0;i<numberRows_;i++) {
3344    CoinExponent expValue=ABC_EXPONENT(region[i]);
3345    if (expValue) {
3346      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
3347        regionIndex[numberNonZero++]=i;
3348      } else {
3349        region[i]=0.0;
3350      }
3351    }
3352  }
3353  regionSparse.setNumElements(numberNonZero);
3354#endif
3355  fromLongArray(static_cast<int>(0));
3356  return numberNonZero;
3357}
3358#if ABC_SMALL<2
3359
3360/* Updates part of column transpose (BTRANU) when densish,
3361   assumes index is sorted i.e. region is correct */
3362void 
3363CoinAbcTypeFactorization::updateColumnTransposeUDensish
3364                        ( CoinIndexedVector * regionSparse,
3365                          CoinSimplexInt smallestIndex) const
3366{
3367  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3368  CoinSimplexInt numberNonZero = regionSparse->getNumElements (  );
3369 
3370  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
3371 
3372  const CoinBigIndex * COIN_RESTRICT startRow = startRowUAddress_;
3373  const CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
3374  const CoinSimplexInt * COIN_RESTRICT numberInRow = numberInRowAddress_;
3375 
3376  const CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
3377 
3378  numberNonZero = 0;
3379  const CoinSimplexInt *  COIN_RESTRICT pivotLinked = pivotLinkedForwardsAddress_;
3380  CoinSimplexInt jRow=smallestIndex ; //pivotLinked[-1];
3381  instrument_start("CoinAbcFactorizationUpdateTransposeUDensish",numberRows_);
3382  while (jRow>=0) {
3383    if (TEST_DOUBLE_NONZERO(region[jRow])) {
3384      CoinFactorizationDouble pivotValue = region[jRow];
3385      if ( !TEST_LESS_THAN_UPDATE_TOLERANCE(pivotValue )  ) {
3386        CoinSimplexInt numberIn = numberInRow[jRow];
3387        instrument_add(numberIn);
3388        if (TEST_INT_NONZERO(numberIn)) {
3389          CoinBigIndex start = startRow[jRow];
3390          CoinBigIndex end = start + numberIn;
3391#ifndef INLINE_IT
3392          for (CoinBigIndex j = start ; j < end; j ++ ) {
3393            CoinSimplexInt iRow = indexColumn[j];
3394            CoinFactorizationDouble value = elementRowU[j];
3395            assert (value);
3396            // i<iRow
3397            region[iRow] -=  value * pivotValue;
3398          }
3399#else
3400          CoinAbcScatterUpdate(end-start,pivotValue,elementRowU+start,indexColumn+start,region);
3401#endif   
3402        } 
3403        regionIndex[numberNonZero++] = jRow;
3404      } else {
3405        region[jRow] = 0.0;
3406      }       
3407    }
3408    jRow=pivotLinked[jRow];
3409  }
3410  //set counts
3411  regionSparse->setNumElements ( numberNonZero );
3412  instrument_end();
3413}
3414/* Updates part of column transpose (BTRANU) when sparse,
3415   assumes index is sorted i.e. region is correct */
3416void 
3417CoinAbcTypeFactorization::updateColumnTransposeUSparse ( 
3418                                                        CoinIndexedVector * regionSparse
3419#if ABC_PARALLEL
3420                                                ,int whichSparse
3421#endif
3422                                                         ) const
3423{
3424  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3425  CoinSimplexInt numberNonZero = regionSparse->getNumElements (  );
3426 
3427  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
3428  const CoinBigIndex * COIN_RESTRICT startRow = startRowUAddress_;
3429  const CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
3430  const CoinSimplexInt * COIN_RESTRICT numberInRow = numberInRowAddress_;
3431 
3432  const CoinFactorizationDouble * COIN_RESTRICT elementRowU = elementRowUAddress_;
3433 
3434  // use sparse_ as temporary area
3435  // mark known to be zero
3436  //printf("PP 0_ %d %s\n",__LINE__,__FILE__);
3437  CoinAbcStack * COIN_RESTRICT stackList = reinterpret_cast<CoinAbcStack *>(sparseAddress_);
3438  CoinSimplexInt * COIN_RESTRICT list = listAddress_;
3439  CoinCheckZero * COIN_RESTRICT mark = markRowAddress_;
3440#if ABC_PARALLEL
3441  //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
3442  if (whichSparse) {
3443    //printf("alternative sparse\n");
3444    int addAmount=whichSparse*sizeSparseArray_;
3445    stackList=reinterpret_cast<CoinAbcStack *>(reinterpret_cast<char *>(stackList)+addAmount);
3446    list=reinterpret_cast<CoinSimplexInt *>(reinterpret_cast<char *>(list)+addAmount);
3447    mark=reinterpret_cast<CoinCheckZero *>(reinterpret_cast<char *>(mark)+addAmount);
3448  }
3449#endif
3450  CoinSimplexInt nList;
3451#if 0 //ndef NDEBUG
3452  for (CoinSimplexInt i=0;i<numberRows_;i++) {
3453#if 0
3454    assert (!mark[i]);
3455#else
3456    if (mark[i]) {
3457      printf("HHHHHHHHHHelp %d %d\n",i,mark[i]);
3458      usleep(10000000);
3459    }
3460#endif
3461  }
3462#endif
3463  nList=0;
3464  for (CoinSimplexInt k=0;k<numberNonZero;k++) {
3465    CoinSimplexInt kPivot=regionIndex[k];
3466    stackList[0].stack=kPivot;
3467    CoinSimplexInt jPivot=kPivot;
3468    CoinBigIndex start=startRow[jPivot];
3469    stackList[0].next=start+numberInRow[jPivot]-1;
3470    stackList[0].start=start;
3471    CoinSimplexInt nStack=1;
3472    while (nStack) {
3473      /* take off stack */
3474      kPivot=stackList[--nStack].stack;
3475      if (mark[kPivot]!=1) {
3476#if 0
3477        CoinBigIndex j=stackList[nStack].next;
3478        CoinBigIndex startThis=stackList[nStack].start;
3479        for (;j>=startThis;j--) {
3480          CoinSimplexInt kPivot3=indexColumn[j--];
3481          if (!mark[kPivot3]) {
3482            kPivot=kPivot3;
3483            break;
3484          }
3485        }
3486        if (j>=startThis) {
3487          /* put back on stack */
3488          stackList[nStack++].next =j;
3489          /* and new one */
3490          CoinBigIndex start=startRow[kPivot];
3491          j=start+numberInRow[kPivot]-1;
3492          stackList[nStack].stack=kPivot;
3493          stackList[nStack].start=start;
3494          mark[kPivot]=2;
3495          stackList[nStack++].next=j;
3496#else
3497        CoinBigIndex j=stackList[nStack].next;
3498        if (j>=stackList[nStack].start) {
3499          CoinSimplexInt kPivot3=indexColumn[j--];
3500          /* put back on stack */
3501          stackList[nStack++].next =j;
3502          if (!mark[kPivot3]) {
3503            /* and new one */
3504            CoinBigIndex start=startRow[kPivot3];
3505            j=start+numberInRow[kPivot3]-1;
3506            stackList[nStack].stack=kPivot3;
3507            stackList[nStack].start=start;
3508            mark[kPivot3]=2;
3509            stackList[nStack++].next=j;
3510          }
3511#endif
3512        } else {
3513          // finished
3514          list[nList++]=kPivot;
3515          mark[kPivot]=1;
3516        }
3517      }
3518    }
3519  }
3520  instrument_start("CoinAbcFactorizationUpdateTransposeUSparse",2*(numberNonZero+nList));
3521  numberNonZero=0;
3522  for (CoinSimplexInt i=nList-1;i>=0;i--) {
3523    CoinSimplexInt iPivot = list[i];
3524    assert (iPivot>=0&&iPivot<numberRows_);
3525    mark[iPivot]=0;
3526    CoinFactorizationDouble pivotValue = region[iPivot];
3527    if ( !TEST_LESS_THAN_UPDATE_TOLERANCE(pivotValue )  ) {
3528      CoinSimplexInt numberIn = numberInRow[iPivot];
3529      instrument_add(numberIn);
3530      if (TEST_INT_NONZERO(numberIn)) {
3531        CoinBigIndex start = startRow[iPivot];
3532#ifndef INLINE_IT
3533        CoinBigIndex end = start + numberIn;
3534        for (CoinBigIndex j=start ; j < end; j ++ ) {
3535          CoinSimplexInt iRow = indexColumn[j];
3536          CoinFactorizationDouble value = elementRowU[j];
3537          region[iRow] -= value * pivotValue;
3538        }
3539#else
3540        CoinAbcScatterUpdate(numberIn,pivotValue,elementRowU+start,indexColumn+start,region);
3541#endif 
3542      }   
3543      regionIndex[numberNonZero++] = iPivot;
3544    } else {
3545      region[iPivot] = 0.0;
3546    }       
3547  }       
3548  //set counts
3549  regionSparse->setNumElements ( numberNonZero );
3550  instrument_end_and_adjust(1.3);
3551}
3552#endif
3553//  updateColumnTransposeU.  Updates part of column transpose (BTRANU)
3554//assumes index is sorted i.e. region is correct
3555//does not sort by sign
3556void
3557CoinAbcTypeFactorization::updateColumnTransposeU ( CoinIndexedVector * regionSparse,
3558                                                   CoinSimplexInt smallestIndex
3559#if ABC_SMALL<2
3560                       , CoinAbcStatistics & statistics
3561#endif
3562#if ABC_PARALLEL
3563                    ,int whichCpu
3564#endif
3565                                                   ) const
3566{
3567#if CILK_CONFLICT>0
3568#if ABC_PARALLEL
3569  // for conflicts
3570#if CILK_CONFLICT>1
3571  printf("file %s line %d which %d\n",__FILE__,__LINE__,whichCpu);
3572#endif
3573  abc_assert((cilk_conflict&(1<<whichCpu))==0);
3574  cilk_conflict |= (1<<whichCpu);
3575#else
3576  abc_assert((cilk_conflict&(1<<0))==0);
3577  cilk_conflict |= (1<<0);
3578#endif
3579#endif
3580#if ABC_SMALL<2
3581  CoinSimplexInt number = regionSparse->getNumElements (  );
3582  if (factorizationStatistics()) {
3583    statistics.numberCounts_++;
3584    statistics.countInput_ += static_cast<CoinSimplexDouble> (number);
3585  }
3586  CoinSimplexInt goSparse;
3587  // Guess at number at end
3588  if (gotUCopy()) {
3589    if (gotSparse()) {
3590      assert (statistics.averageAfterU_);
3591      CoinSimplexInt newNumber = static_cast<CoinSimplexInt> (number*statistics.averageAfterU_*twiddleBtranFactor1());
3592      if (newNumber< sparseThreshold_)
3593        goSparse = 2;
3594      else
3595        goSparse = 0;
3596    } else {
3597      goSparse=0;
3598    }
3599#if ABC_SMALL>=0
3600  } else {
3601    goSparse=-1;
3602#endif
3603  }
3604  //CoinIndexedVector temp(*regionSparse);
3605  switch (goSparse) {
3606#if ABC_SMALL>=0
3607  case -1: // no row copy
3608    updateColumnTransposeUByColumn(regionSparse,smallestIndex);
3609    break;
3610#endif
3611  case 0: // densish
3612    updateColumnTransposeUDensish(regionSparse,smallestIndex);
3613    break;
3614  case 2: // sparse
3615    {
3616      updateColumnTransposeUSparse(regionSparse
3617#if ABC_PARALLEL
3618                                      ,whichCpu
3619#endif
3620                                   );
3621    }
3622    break;
3623  }
3624  if (factorizationStatistics()) 
3625    statistics.countAfterU_ += static_cast<CoinSimplexDouble> (regionSparse->getNumElements());
3626#else
3627  updateColumnTransposeUByColumn(regionSparse,smallestIndex);
3628#endif
3629#if CILK_CONFLICT>0
3630#if ABC_PARALLEL
3631  // for conflicts
3632  abc_assert((cilk_conflict&(1<<whichCpu))!=0);
3633  cilk_conflict &= ~(1<<whichCpu);
3634#else
3635  abc_assert((cilk_conflict&(1<<0))!=0);
3636  cilk_conflict &= ~(1<<0);
3637#endif
3638#endif
3639}
3640
3641/*  updateColumnTransposeLDensish. 
3642    Updates part of column transpose (BTRANL) dense by column */
3643void
3644CoinAbcTypeFactorization::updateColumnTransposeLDensish
3645     ( CoinIndexedVector * regionSparse ) const
3646{
3647  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3648  CoinSimplexInt base;
3649  CoinSimplexInt first = -1;
3650  const CoinSimplexInt * COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
3651  //const CoinSimplexInt * COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
3652 
3653#if ABC_SMALL<3
3654  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
3655  CoinSimplexInt numberNonZero=0;
3656#endif
3657  //scan
3658  for (first=numberRows_-1;first>=0;first--) {
3659#ifndef ABC_ORDERED_FACTORIZATION
3660    if (region[pivotLOrder[first]]) 
3661      break;
3662#else
3663    if (region[first]) 
3664      break;
3665#endif
3666  }
3667  instrument_start("CoinAbcFactorizationUpdateTransposeLDensish",first);
3668  if ( first >= 0 ) {
3669    base = baseL_;
3670    const CoinBigIndex * COIN_RESTRICT startColumn = startColumnLAddress_;
3671    const CoinSimplexInt * COIN_RESTRICT indexRow = indexRowLAddress_;
3672    const CoinFactorizationDouble * COIN_RESTRICT element = elementLAddress_;
3673    CoinSimplexInt last = baseL_ + numberL_;
3674   
3675    if ( first >= last ) {
3676      first = last - 1;
3677    } 
3678    CoinBigIndex end = startColumn[first+1];
3679    for (CoinSimplexInt k = first ; k >= base; k-- ) {
3680#ifndef ABC_ORDERED_FACTORIZATION
3681      CoinSimplexInt i=pivotLOrder[k];
3682#else
3683      CoinSimplexInt i=k;
3684#endif
3685      CoinFactorizationDouble pivotValue = region[i];
3686#ifndef INLINE_IT2
3687      CoinBigIndex j=end-1;
3688      end=startColumn[k];
3689      instrument_add(j-end);
3690      for (  ; j >= end; j-- ) {
3691        CoinSimplexInt iRow = indexRow[j];
3692        CoinFactorizationDouble value = element[j];
3693        pivotValue -= value * region[iRow];
3694      }       
3695#else
3696      CoinBigIndex start=startColumn[k];
3697      instrument_add(end-start);
3698      pivotValue+=CoinAbcGatherUpdate(end-start,element+start,indexRow+start,region);
3699      end=start;
3700#endif
3701#if ABC_SMALL<3
3702      if ( !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue )  ) {
3703        region[i] = pivotValue;
3704        regionIndex[numberNonZero++] = i;
3705      } else { 
3706        region[i] = 0.0;
3707      }
3708#else
3709        region[i] = pivotValue;
3710#endif       
3711    }       
3712#if ABC_SMALL<3
3713    //may have stopped early
3714    if ( first < base ) {
3715      base = first + 1;
3716    }
3717   
3718    for (CoinSimplexInt k = base-1 ; k >= 0; k-- ) {
3719#ifndef ABC_ORDERED_FACTORIZATION
3720      CoinSimplexInt i=pivotLOrder[k];
3721#else
3722      CoinSimplexInt i=k;
3723#endif
3724      CoinExponent expValue=ABC_EXPONENT(region[i]);
3725      if (expValue) {
3726        if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
3727          regionIndex[numberNonZero++] = i;
3728        } else {
3729          region[i] = 0.0;
3730        }       
3731      }
3732    }
3733#endif
3734  } 
3735#if ABC_SMALL<3
3736  //set counts
3737  regionSparse->setNumElements ( numberNonZero );
3738#endif
3739  instrument_end();
3740}
3741#if ABC_SMALL<2
3742/*  updateColumnTransposeLByRow.
3743    Updates part of column transpose (BTRANL) densish but by row */
3744void
3745CoinAbcTypeFactorization::updateColumnTransposeLByRow
3746    ( CoinIndexedVector * regionSparse ) const
3747{
3748  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3749  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
3750  CoinSimplexInt numberNonZero;
3751 
3752  // use row copy of L
3753  const CoinFactorizationDouble *  COIN_RESTRICT element = elementByRowLAddress_;
3754  const CoinBigIndex *  COIN_RESTRICT startRow = startRowLAddress_;
3755  const CoinSimplexInt *  COIN_RESTRICT column = indexColumnLAddress_;
3756  const CoinSimplexInt * COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
3757  numberNonZero=0;
3758  instrument_start("CoinAbcFactorizationUpdateTransposeLByRow",numberRows_+numberSlacks_);
3759  // down to slacks
3760  for (CoinSimplexInt k = numberRows_-1 ; k >= numberSlacks_; k-- ) {
3761#ifndef ABC_ORDERED_FACTORIZATION
3762    CoinSimplexInt i=pivotLOrder[k];
3763#else
3764    CoinSimplexInt i=k;
3765#endif
3766    if (TEST_DOUBLE_NONZERO(region[i])) {
3767      CoinFactorizationDouble pivotValue = region[i];
3768      if ( !TEST_LESS_THAN_TOLERANCE(pivotValue )  ) {
3769        regionIndex[numberNonZero++] = i;
3770        CoinBigIndex start=startRow[i];
3771        CoinBigIndex end=startRow[i+1];
3772        instrument_add(end-start);
3773        if (TEST_INT_NONZERO(end-start)) {
3774#ifndef INLINE_IT
3775          for (CoinBigIndex j = end-1;j >= start; j--) {
3776            CoinSimplexInt iRow = column[j];
3777            CoinFactorizationDouble value = element[j];
3778            region[iRow] -= pivotValue*value;
3779          }
3780#else
3781          CoinAbcScatterUpdate(end-start,pivotValue,element+start,column+start,region);
3782#endif
3783        }
3784      } else {
3785        region[i] = 0.0;
3786      }     
3787    }
3788  }
3789  for (CoinSimplexInt k = numberSlacks_-1 ; k >= 0; k-- ) {
3790#ifndef ABC_ORDERED_FACTORIZATION
3791    CoinSimplexInt i=pivotLOrder[k];
3792#else
3793    CoinSimplexInt i=k;
3794#endif
3795    CoinExponent expValue=ABC_EXPONENT(region[i]);
3796    if (expValue) {
3797      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
3798        regionIndex[numberNonZero++]=i;
3799      } else {
3800        region[i] = 0.0;
3801      }     
3802    }
3803  }
3804  //set counts
3805  regionSparse->setNumElements ( numberNonZero );
3806  instrument_end();
3807}
3808/*  updateColumnTransposeLSparse.
3809    Updates part of column transpose (BTRANL) sparse */
3810void
3811CoinAbcTypeFactorization::updateColumnTransposeLSparse ( CoinIndexedVector * regionSparse
3812#if ABC_PARALLEL
3813                                                ,int whichSparse
3814#endif
3815                                                         ) const
3816{
3817  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3818  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
3819  CoinSimplexInt numberNonZero = regionSparse->getNumElements (  );
3820 
3821  // use row copy of L
3822  const CoinFactorizationDouble *  COIN_RESTRICT element = elementByRowLAddress_;
3823  const CoinBigIndex *  COIN_RESTRICT startRow = startRowLAddress_;
3824  const CoinSimplexInt *  COIN_RESTRICT column = indexColumnLAddress_;
3825  // use sparse_ as temporary area
3826  // mark known to be zero
3827  //printf("PP 0_ %d %s\n",__LINE__,__FILE__);
3828  CoinAbcStack * COIN_RESTRICT stackList = reinterpret_cast<CoinAbcStack *>(sparseAddress_);
3829  CoinSimplexInt * COIN_RESTRICT list = listAddress_;
3830  CoinCheckZero * COIN_RESTRICT mark = markRowAddress_;
3831#if ABC_PARALLEL
3832  //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
3833  if (whichSparse) {
3834    int addAmount=whichSparse*sizeSparseArray_;
3835    //printf("alternative sparse\n");
3836    stackList=reinterpret_cast<CoinAbcStack *>(reinterpret_cast<char *>(stackList)+addAmount);
3837    list=reinterpret_cast<CoinSimplexInt *>(reinterpret_cast<char *>(list)+addAmount);
3838    mark=reinterpret_cast<CoinCheckZero *>(reinterpret_cast<char *>(mark)+addAmount);
3839  }
3840#endif
3841#ifndef NDEBUG
3842  for (CoinSimplexInt i=0;i<numberRows_;i++) {
3843#if 0
3844    assert (!mark[i]);
3845#else
3846    if (mark[i]) {
3847      printf("HHHHHHHHHHelp2 %d %d\n",i,mark[i]);
3848      usleep(10000000);
3849    }
3850#endif
3851  }
3852#endif
3853  CoinSimplexInt nList;
3854  CoinSimplexInt number = numberNonZero;
3855  nList=0;
3856  for (CoinSimplexInt k=0;k<number;k++) {
3857    CoinSimplexInt kPivot=regionIndex[k];
3858    //CoinSimplexInt iPivot=pivotLBackwardOrder[kPivot];
3859    if(!mark[kPivot]&&region[kPivot]) {
3860      stackList[0].stack=kPivot;
3861      stackList[0].start=startRow[kPivot];
3862      CoinBigIndex j=startRow[kPivot+1]-1;
3863      CoinSimplexInt nStack=0;
3864      while (nStack>=0) {
3865        /* take off stack */
3866        if (j>=stackList[nStack].start) {
3867          CoinSimplexInt jPivot=column[j--];
3868          /* put back on stack */
3869          stackList[nStack].next =j;
3870          if (!mark[jPivot]) {
3871            /* and new one */
3872            kPivot=jPivot;
3873            //iPivot=pivotLBackwardOrder[kPivot];
3874            j = startRow[kPivot+1]-1;
3875            stackList[++nStack].stack=kPivot;
3876            mark[kPivot]=1;
3877            stackList[nStack].start=startRow[kPivot];
3878            stackList[nStack].next=j;
3879          }
3880        } else {
3881          /* finished so mark */
3882          list[nList++]=kPivot;
3883          mark[kPivot]=1;
3884          --nStack;
3885          if (nStack>=0) {
3886            kPivot=stackList[nStack].stack;
3887            j=stackList[nStack].next;
3888          }
3889        }
3890      }
3891    }
3892  }
3893  instrument_start("CoinAbcFactorizationUpdateTransposeLSparse",2*(numberNonZero+nList));
3894  numberNonZero=0;
3895  for (CoinSimplexInt i=nList-1;i>=0;i--) {
3896    CoinSimplexInt iPivot = list[i];
3897    //CoinSimplexInt kPivot=pivotLOrder[iPivot];
3898    mark[iPivot]=0;
3899    CoinFactorizationDouble pivotValue = region[iPivot];
3900    if ( !TEST_LESS_THAN_TOLERANCE(pivotValue )  ) {
3901      regionIndex[numberNonZero++] = iPivot;
3902      CoinBigIndex start=startRow[iPivot];
3903      CoinBigIndex end=startRow[iPivot+1];
3904      instrument_add(end-start);
3905#ifndef INLINE_IT
3906      CoinBigIndex j;
3907      for ( j = start; j < end; j ++ ) {
3908        CoinSimplexInt iRow = column[j];
3909        CoinFactorizationDouble value = element[j];
3910        region[iRow] -= value * pivotValue;
3911      }
3912#else
3913      CoinAbcScatterUpdate(end-start,pivotValue,element+start,column+start,region);
3914#endif
3915    } else {
3916      region[iPivot]=0.0;
3917    }
3918  }
3919  //set counts
3920  regionSparse->setNumElements ( numberNonZero );
3921  instrument_end_and_adjust(1.3);
3922}
3923#endif
3924//  updateColumnTransposeL.  Updates part of column transpose (BTRANL)
3925void
3926CoinAbcTypeFactorization::updateColumnTransposeL ( CoinIndexedVector * regionSparse
3927#if ABC_SMALL<2
3928                       , CoinAbcStatistics & statistics
3929#endif
3930#if ABC_PARALLEL
3931                                      ,int whichSparse
3932#endif
3933                                                   ) const
3934{
3935#if CILK_CONFLICT>0
3936#if ABC_PARALLEL
3937  // for conflicts
3938#if CILK_CONFLICT>1
3939  printf("file %s line %d which %d\n",__FILE__,__LINE__,whichSparse);
3940#endif
3941  abc_assert((cilk_conflict&(1<<whichSparse))==0);
3942  cilk_conflict |= (1<<whichSparse);
3943#else
3944  abc_assert((cilk_conflict&(1<<0))==0);
3945  cilk_conflict |= (1<<0);
3946#endif
3947#endif
3948#if ABC_SMALL<3
3949  CoinSimplexInt number = regionSparse->getNumElements (  );
3950#endif
3951  if (!numberL_
3952#if ABC_SMALL<4
3953      &&!numberDense_
3954#endif
3955      ) {
3956#if ABC_SMALL<3
3957    if (number>=numberRows_) {
3958      // need scan
3959      regionSparse->setNumElements(0);
3960      scan(regionSparse);
3961    }
3962#endif
3963#if CILK_CONFLICT>0
3964#if ABC_PARALLEL
3965    // for conflicts
3966    abc_assert((cilk_conflict&(1<<whichSparse))!=0);
3967    cilk_conflict &= ~(1<<whichSparse);
3968#else
3969    abc_assert((cilk_conflict&(1<<0))!=0);
3970    cilk_conflict &= ~(1<<0);
3971#endif
3972#endif
3973    return;
3974  }
3975#if ABC_SMALL<2
3976  CoinSimplexInt goSparse;
3977  // Guess at number at end
3978  // we may need to rethink on dense
3979  if (gotLCopy()) {
3980    assert (statistics.averageAfterL_);
3981    CoinSimplexInt newNumber = static_cast<CoinSimplexInt> (number*statistics.averageAfterL_*twiddleBtranFactor1());
3982    if (newNumber< sparseThreshold_)
3983      goSparse = 2;
3984    else if (2*newNumber< numberRows_)
3985      goSparse = 0;
3986    else
3987      goSparse = -1;
3988  } else {
3989    goSparse=-1;
3990  }
3991#endif
3992#if ABC_SMALL<4
3993#if ABC_DENSE_CODE>0
3994  if (numberDense_) {
3995    //take off list
3996    CoinSimplexInt lastSparse = numberRows_-numberDense_;
3997    CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
3998    const CoinSimplexInt * COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
3999    CoinFactorizationDouble *  COIN_RESTRICT denseArea = denseAreaAddress_;
4000    CoinFactorizationDouble * COIN_RESTRICT denseRegion = 
4001      denseArea+leadingDimension_*numberDense_;
4002    CoinSimplexInt *  COIN_RESTRICT densePermute=
4003      reinterpret_cast<CoinSimplexInt *>(denseRegion+FACTOR_CPU*numberDense_);
4004#if ABC_PARALLEL
4005    //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
4006    if (whichSparse)
4007      denseRegion+=whichSparse*numberDense_;
4008#endif
4009    CoinFactorizationDouble * COIN_RESTRICT densePut = 
4010      denseRegion-lastSparse;
4011    CoinZeroN(denseRegion,numberDense_);
4012    bool doDense=false;
4013#if ABC_SMALL<3
4014    CoinSimplexInt number = regionSparse->getNumElements();
4015    CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
4016    const CoinSimplexInt * COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
4017    CoinSimplexInt i=0;
4018#if ABC_SMALL<0
4019    assert (number<=numberRows_);
4020#else
4021    if (number<=numberRows_) {
4022#endif
4023      while (i<number) {
4024        CoinSimplexInt iRow = regionIndex[i];
4025#ifndef ABC_ORDERED_FACTORIZATION
4026        CoinSimplexInt jRow = pivotLBackwardOrder[iRow];
4027#else
4028        CoinSimplexInt jRow=iRow;
4029#endif
4030        if (jRow>=lastSparse) {
4031          doDense=true;
4032          regionIndex[i] = regionIndex[--number];
4033          densePut[jRow]=region[iRow];
4034          region[iRow]=0.0;
4035        } else {
4036          i++;
4037        }
4038      }
4039#endif
4040#if ABC_SMALL>=0
4041#if ABC_SMALL<3
4042    } else {
4043#endif
4044      for (CoinSimplexInt jRow=lastSparse;jRow<numberRows_;jRow++) {
4045        CoinSimplexInt iRow = pivotLOrder[jRow];
4046        if (region[iRow]) {
4047          doDense=true;
4048          densePut[jRow]=region[iRow];
4049          region[iRow]=0.0;
4050        }
4051      }
4052#if ABC_SMALL<3
4053      // numbers are all wrong - do a scan
4054      regionSparse->setNumElements(0);
4055      scan(regionSparse);
4056      number=regionSparse->getNumElements();
4057    }
4058#endif
4059#endif
4060    if (doDense) {
4061      instrument_do("CoinAbcFactorizationDenseTranspose",0.25*numberDense_*numberDense_);
4062#if ABC_DENSE_CODE==1
4063#ifndef CLAPACK
4064      char trans = 'T';
4065      CoinSimplexInt ione=1;
4066      CoinSimplexInt info;
4067      F77_FUNC(dgetrs,DGETRS)(&trans,&numberDense_,&ione,denseArea,&leadingDimension_,
4068                              densePermute,denseRegion,&numberDense_,&info,1);
4069#else
4070      clapack_dgetrs(CblasColMajor, CblasTrans,numberDense_,1,
4071                     denseArea, leadingDimension_,densePermute,denseRegion,numberDense_);
4072#endif
4073#elif ABC_DENSE_CODE==2
4074      CoinAbcDgetrs('T',numberDense_,denseArea,denseRegion);
4075      short * COIN_RESTRICT forBtran = reinterpret_cast<short *>(densePermute+numberDense_)+numberDense_-lastSparse;
4076      pivotLOrder += lastSparse; // adjust
4077#endif
4078      for (CoinSimplexInt i=lastSparse;i<numberRows_;i++) {
4079        CoinSimplexDouble value = densePut[i];
4080        if (value) {
4081          if (!TEST_LESS_THAN_TOLERANCE(value)) {
4082#if ABC_DENSE_CODE!=2
4083            CoinSimplexInt iRow=pivotLOrder[i];
4084#else
4085            CoinSimplexInt iRow=pivotLOrder[forBtran[i]];
4086#endif
4087            region[iRow]=value;
4088#if ABC_SMALL<3
4089            regionIndex[number++] = iRow;
4090#endif
4091          }
4092        }
4093      }
4094#if ABC_SMALL<3
4095      regionSparse->setNumElements(number);
4096#endif
4097    } 
4098  } 
4099  //printRegion(*regionSparse,"After BtranL");
4100#endif
4101#endif
4102#if ABC_SMALL<2
4103  switch (goSparse) {
4104  case -1: // No row copy
4105    updateColumnTransposeLDensish(regionSparse);
4106    break;
4107  case 0: // densish but by row
4108    updateColumnTransposeLByRow(regionSparse);
4109    break;
4110  case 2: // sparse
4111    updateColumnTransposeLSparse(regionSparse
4112#if ABC_PARALLEL
4113                                      ,whichSparse
4114#endif
4115                                 );
4116    break;
4117  }
4118#else
4119  updateColumnTransposeLDensish(regionSparse);
4120#endif
4121#if CILK_CONFLICT>0
4122#if ABC_PARALLEL
4123  // for conflicts
4124  abc_assert((cilk_conflict&(1<<whichSparse))!=0);
4125  cilk_conflict &= ~(1<<whichSparse);
4126#else
4127  abc_assert((cilk_conflict&(1<<0))!=0);
4128  cilk_conflict &= ~(1<<0);
4129#endif
4130#endif
4131}
4132#if ABC_SMALL>=0
4133// Updates part of column transpose (BTRANR) when dense
4134void 
4135CoinAbcTypeFactorization::updateColumnTransposeRDensish( CoinIndexedVector * regionSparse ) const
4136{
4137  const CoinBigIndex *  COIN_RESTRICT startColumn = startColumnRAddress_-numberRows_;
4138  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
4139#if ABC_SMALL<3
4140#ifdef ABC_DEBUG
4141  regionSparse->checkClean();
4142#endif
4143#endif
4144  instrument_start("CoinAbcFactorizationUpdateTransposeRDensish",numberRows_);
4145  CoinSimplexInt last = numberRowsExtra_-1;
4146  const CoinSimplexInt * COIN_RESTRICT indexRow = indexRowRAddress_;
4147  const CoinFactorizationDouble * COIN_RESTRICT element = elementRAddress_;
4148  const CoinSimplexInt *  COIN_RESTRICT pivotRowBack = pivotColumnAddress_;
4149  for (CoinSimplexInt i = last ; i >= numberRows_; i-- ) {
4150    CoinSimplexInt putRow=pivotRowBack[i];
4151    CoinFactorizationDouble pivotValue = region[putRow];
4152    if ( pivotValue ) {
4153      CoinBigIndex start=startColumn[i];
4154      CoinBigIndex end=startColumn[i+1];
4155      instrument_add(end-start);
4156#ifndef INLINE_IT
4157      for (CoinBigIndex j = start; j < end; j++ ) {
4158        CoinFactorizationDouble value = element[j];
4159        CoinSimplexInt iRow = indexRow[j];
4160        assert (iRow<numberRows_);
4161        region[iRow] -= value * pivotValue;
4162      }
4163#else
4164      CoinAbcScatterUpdate(end-start,pivotValue,element+start,indexRow+start,region);
4165#endif
4166      region[putRow] = pivotValue;
4167    }
4168  }
4169  instrument_end();
4170}
4171#endif
4172#if ABC_SMALL<2
4173// Updates part of column transpose (BTRANR) when sparse
4174void 
4175CoinAbcTypeFactorization::updateColumnTransposeRSparse
4176( CoinIndexedVector * regionSparse ) const
4177{
4178  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(regionSparse);
4179  CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
4180  CoinSimplexInt numberNonZero = regionSparse->getNumElements (  );
4181  instrument_start("CoinAbcFactorizationUpdateTransposeRSparse",numberRows_);
4182
4183#if ABC_SMALL<3
4184#ifdef ABC_DEBUG
4185  regionSparse->checkClean();
4186#endif
4187#endif
4188  CoinSimplexInt last = numberRowsExtra_-1;
4189 
4190  const CoinSimplexInt * COIN_RESTRICT indexRow = indexRowRAddress_;
4191  const CoinFactorizationDouble * COIN_RESTRICT element = elementRAddress_;
4192  const CoinBigIndex *  COIN_RESTRICT startColumn = startColumnRAddress_-numberRows_;
4193  //move using lookup
4194  const CoinSimplexInt *  COIN_RESTRICT pivotRowBack = pivotColumnAddress_;
4195  // still need to do in correct order (for now)
4196  for (CoinSimplexInt i = last ; i >= numberRows_; i-- ) {
4197    CoinSimplexInt putRow = pivotRowBack[i];
4198    CoinFactorizationDouble pivotValue = region[putRow];
4199    if ( pivotValue ) {
4200      instrument_add(startColumn[i+1]-startColumn[i]);
4201      for (CoinBigIndex j = startColumn[i]; j < startColumn[i+1]; j++ ) {
4202        CoinFactorizationDouble value = element[j];
4203        CoinSimplexInt iRow = indexRow[j];
4204        assert (iRow<numberRows_);
4205        bool oldValue =  (TEST_DOUBLE_REALLY_NONZERO(region[iRow]));
4206        CoinFactorizationDouble newValue = region[iRow] - value * pivotValue;
4207        if (oldValue) {
4208          if (newValue) 
4209            region[iRow]=newValue;
4210          else
4211            region[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
4212        } else if (!TEST_LESS_THAN_TOLERANCE_REGISTER(newValue)) {
4213          region[iRow] = newValue;
4214          regionIndex[numberNonZero++]=iRow;
4215        }
4216      }
4217      region[putRow] = pivotValue;
4218    }
4219  }
4220  instrument_end();
4221  regionSparse->setNumElements(numberNonZero);
4222}
4223#endif
4224//  updateColumnTransposeR.  Updates part of column (FTRANR)
4225void
4226CoinAbcTypeFactorization::updateColumnTransposeR ( CoinIndexedVector * regionSparse
4227#if ABC_SMALL<2
4228                                                   , CoinAbcStatistics & statistics
4229#endif
4230                                                   ) const
4231{
4232  if (numberRowsExtra_==numberRows_)
4233    return;
4234#if ABC_SMALL<3
4235  CoinSimplexInt numberNonZero = regionSparse->getNumElements (  );
4236  if (numberNonZero) {
4237#endif
4238#if ABC_SMALL<3
4239    const CoinBigIndex *  COIN_RESTRICT startColumn = startColumnRAddress_-numberRows_;
4240#endif
4241    // Size of R
4242    instrument_do("CoinAbcFactorizationTransposeR",startColumnRAddress_[numberR_]);
4243#if ABC_SMALL<2
4244    if (startColumn[numberRowsExtra_]) {
4245#if ABC_SMALL>=0
4246      if (numberNonZero < ((sparseThreshold_<<2)*twiddleBtranFactor2())
4247          ||(!numberL_&&gotLCopy()&&gotRCopy())) {
4248#endif
4249        updateColumnTransposeRSparse ( regionSparse );
4250        if (factorizationStatistics()) 
4251          statistics.countAfterR_ += regionSparse->getNumElements();
4252#if ABC_SMALL>=0
4253      } else {
4254        updateColumnTransposeRDensish ( regionSparse );
4255        // we have lost indices
4256        // make sure won't try and go sparse again
4257        if (factorizationStatistics()) 
4258          statistics.countAfterR_ += CoinMin((numberNonZero<<1),numberRows_);
4259        // temp +2 (was +1)
4260        regionSparse->setNumElements (numberRows_+2);
4261      }
4262#endif
4263    } else {
4264      if (factorizationStatistics()) 
4265        statistics.countAfterR_ += numberNonZero;
4266    }
4267#else
4268    updateColumnTransposeRDensish ( regionSparse );
4269#if ABC_SMALL<3
4270    // we have lost indices
4271    // make sure won't try and go sparse again
4272    // temp +2 (was +1)
4273    regionSparse->setNumElements (numberRows_+2);
4274#endif
4275#endif
4276#if ABC_SMALL<3
4277  }
4278#endif
4279}
4280// Update partial Ftran by R update
4281void 
4282CoinAbcTypeFactorization::updatePartialUpdate(CoinIndexedVector & partialUpdate)
4283{
4284  CoinSimplexDouble * COIN_RESTRICT region = partialUpdate.denseVector (  );
4285  CoinSimplexInt * COIN_RESTRICT regionIndex = partialUpdate.getIndices (  );
4286  CoinSimplexInt numberNonZero = partialUpdate.getNumElements (  );
4287  const CoinSimplexInt *  COIN_RESTRICT indexRow = indexRowRAddress_;
4288  const CoinFactorizationDouble *  COIN_RESTRICT element = elementRAddress_;
4289  const CoinSimplexInt * pivotRowBack = pivotColumnAddress_;
4290  const CoinBigIndex *  COIN_RESTRICT startColumn = startColumnRAddress_-numberRows_;
4291  CoinBigIndex start = startColumn[numberRowsExtra_-1];
4292  CoinBigIndex end = startColumn[numberRowsExtra_];
4293  CoinSimplexInt iRow =pivotRowBack[numberRowsExtra_-1];
4294  CoinFactorizationDouble pivotValue = region[iRow];
4295  for (CoinBigIndex j = start; j < end; j ++ ) {
4296    CoinFactorizationDouble value = element[j];
4297    CoinSimplexInt jRow = indexRow[j];
4298    value *= region[jRow];
4299    pivotValue -= value;
4300  }
4301  if ( !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue ) ) {
4302    if (!region[iRow])
4303      regionIndex[numberNonZero++] = iRow;
4304    region[iRow] = pivotValue;
4305  } else {
4306    if (region[iRow])
4307      region[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
4308  }
4309  partialUpdate.setNumElements(numberNonZero);
4310}
4311#if FACTORIZATION_STATISTICS
4312#ifdef ABC_JUST_ONE_FACTORIZATION
4313#define FACTORS_HERE
4314#endif
4315#ifdef FACTORS_HERE
4316double ftranTwiddleFactor1X=1.0;
4317double ftranTwiddleFactor2X=1.0;
4318double ftranFTTwiddleFactor1X=1.0;
4319double ftranFTTwiddleFactor2X=1.0;
4320double btranTwiddleFactor1X=1.0;
4321double btranTwiddleFactor2X=1.0;
4322double ftranFullTwiddleFactor1X=1.0;
4323double ftranFullTwiddleFactor2X=1.0;
4324double btranFullTwiddleFactor1X=1.0;
4325double btranFullTwiddleFactor2X=1.0;
4326double denseThresholdX=31;
4327double twoThresholdX=1000;
4328double minRowsSparse=300;
4329double largeRowsSparse=10000;
4330double mediumRowsDivider=6;
4331double mediumRowsMinCount=500;
4332double largeRowsCount=1000;
4333#else
4334extern double ftranTwiddleFactor1X;
4335extern double ftranTwiddleFactor2X;
4336extern double ftranFTTwiddleFactor1X;
4337extern double ftranFTTwiddleFactor2X;
4338extern double btranTwiddleFactor1X;
4339extern double btranTwiddleFactor2X;
4340extern double ftranFullTwiddleFactor1X;
4341extern double ftranFullTwiddleFactor2X;
4342extern double btranFullTwiddleFactor1X;
4343extern double btranFullTwiddleFactor2X;
4344extern double denseThresholdX;
4345extern double twoThresholdX;
4346extern double minRowsSparse;
4347extern double largeRowsSparse;
4348extern double mediumRowsDivider;
4349extern double mediumRowsMinCount;
4350extern double largeRowsCount;
4351#endif
4352#else
4353#define minRowsSparse 300
4354#define largeRowsSparse 10000
4355#define mediumRowsDivider 6
4356#define mediumRowsMinCount 500
4357#define largeRowsCount 1000
4358#endif
4359//  makes a row copy of L
4360void
4361CoinAbcTypeFactorization::goSparse2 ( )
4362{
4363#if ABC_SMALL<2
4364#if ABC_SMALL>=0
4365  if (numberRows_>minRowsSparse) {
4366#endif
4367    if (numberRows_<largeRowsSparse) {
4368      sparseThreshold_=CoinMin(numberRows_/mediumRowsDivider,mediumRowsMinCount);
4369    } else {
4370      sparseThreshold_=CoinMax(largeRowsCount,numberRows_>>3);
4371    }
4372#if FACTORIZATION_STATISTICS
4373    ftranTwiddleFactor1_=ftranTwiddleFactor1X;
4374    ftranTwiddleFactor2_=ftranTwiddleFactor2X;
4375    ftranFTTwiddleFactor1_=ftranFTTwiddleFactor1X;
4376    ftranFTTwiddleFactor2_=ftranFTTwiddleFactor2X;
4377    btranTwiddleFactor1_=btranTwiddleFactor1X;
4378    btranTwiddleFactor2_=btranTwiddleFactor2X;
4379    ftranFullTwiddleFactor1_=ftranFullTwiddleFactor1X;
4380    ftranFullTwiddleFactor2_=ftranFullTwiddleFactor2X;
4381    btranFullTwiddleFactor1_=btranFullTwiddleFactor1X;
4382    btranFullTwiddleFactor2_=btranFullTwiddleFactor2X;
4383#endif
4384    setYesGotLCopy();
4385    setYesGotUCopy();
4386    setYesGotSparse();
4387    // allow for stack, list, next, starts and char map of mark
4388    CoinSimplexInt nRowIndex = (numberRows_+CoinSizeofAsInt(CoinSimplexInt)-1)/
4389      CoinSizeofAsInt(char);
4390    CoinSimplexInt nInBig = static_cast<CoinSimplexInt>(sizeof(CoinBigIndex)/sizeof(CoinSimplexInt));
4391    assert (nInBig>=1);
4392    sizeSparseArray_ = (2+2*nInBig)*numberRows_ + nRowIndex ;
4393    sparse_.conditionalNew(FACTOR_CPU*sizeSparseArray_);
4394    sparseAddress_ = sparse_.array();
4395    // zero out mark
4396    CoinAbcStack * stackList = reinterpret_cast<CoinAbcStack *>(sparseAddress_);
4397    listAddress_ = reinterpret_cast<CoinSimplexInt *>(stackList+numberRows_);
4398    markRowAddress_ = reinterpret_cast<CoinCheckZero *> (listAddress_ + numberRows_);
4399    assert(reinterpret_cast<CoinCheckZero *>(sparseAddress_+(2+2*nInBig)*numberRows_)==markRowAddress_);
4400    CoinAbcMemset0(markRowAddress_,nRowIndex);
4401#if ABC_PARALLEL
4402    //printf("PP 0__ %d %s\n",__LINE__,__FILE__);
4403    // convert to bytes
4404    sizeSparseArray_*=sizeof(CoinSimplexInt);
4405    char * mark=reinterpret_cast<char *>(reinterpret_cast<char *>(markRowAddress_)+sizeSparseArray_);
4406    for (int i=1;i<FACTOR_CPU;i++) {
4407      CoinAbcMemset0(mark,nRowIndex);
4408      mark=reinterpret_cast<char *>(reinterpret_cast<char *>(mark)+sizeSparseArray_);
4409    }
4410#endif
4411    elementByRowL_.conditionalDelete();
4412    indexColumnL_.conditionalDelete();
4413    startRowL_.conditionalNew(numberRows_+1);
4414    if (lengthAreaL_) {
4415      elementByRowL_.conditionalNew(lengthAreaL_);
4416      indexColumnL_.conditionalNew(lengthAreaL_);
4417    }
4418    elementByRowLAddress_=elementByRowL_.array();
4419    indexColumnLAddress_=indexColumnL_.array();
4420    startRowLAddress_=startRowL_.array();
4421    // counts
4422    CoinBigIndex * COIN_RESTRICT startRowL = startRowLAddress_;
4423    CoinZeroN(startRowL,numberRows_);
4424    const CoinBigIndex * startColumnL = startColumnLAddress_;
4425    CoinFactorizationDouble * COIN_RESTRICT elementL = elementLAddress_;
4426    const CoinSimplexInt * indexRowL = indexRowLAddress_;
4427    for (CoinSimplexInt i=baseL_;i<baseL_+numberL_;i++) {
4428      for (CoinBigIndex j=startColumnL[i];j<startColumnL[i+1];j++) {
4429        CoinSimplexInt iRow = indexRowL[j];
4430        startRowL[iRow]++;
4431      }
4432    }
4433    const CoinSimplexInt * COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
4434    // convert count to lasts
4435    CoinBigIndex count=0;
4436    for (CoinSimplexInt i=0;i<numberRows_;i++) {
4437      CoinSimplexInt numberInRow=startRowL[i];
4438      count += numberInRow;
4439      startRowL[i]=count;
4440    }
4441    startRowL[numberRows_]=count;
4442    // now insert
4443    CoinFactorizationDouble * COIN_RESTRICT elementByRowL = elementByRowLAddress_;
4444    CoinSimplexInt * COIN_RESTRICT indexColumnL = indexColumnLAddress_;
4445    for (CoinSimplexInt i=baseL_+numberL_-1;i>=baseL_;i--) {
4446#ifndef ABC_ORDERED_FACTORIZATION
4447      CoinSimplexInt iPivot=pivotLOrder[i];
4448#else
4449      CoinSimplexInt iPivot=i;
4450#endif
4451      for (CoinBigIndex j=startColumnL[i];j<startColumnL[i+1];j++) {
4452        CoinSimplexInt iRow = indexRowL[j];
4453        CoinBigIndex start = startRowL[iRow]-1;
4454        startRowL[iRow]=start;
4455        elementByRowL[start]=elementL[j];
4456        indexColumnL[start]=iPivot;
4457      }
4458    }
4459#if ABC_SMALL>=0
4460  } else {
4461    sparseThreshold_=0;
4462    setNoGotLCopy();
4463    //setYesGotUCopy();
4464    setNoGotUCopy();
4465    setNoGotSparse();
4466  }
4467#endif
4468#endif
4469}
4470
4471//  set sparse threshold
4472void
4473CoinAbcTypeFactorization::sparseThreshold ( CoinSimplexInt /*value*/ ) 
4474{
4475  return;
4476#if 0
4477  if (value>0&&sparseThreshold_) {
4478    sparseThreshold_=value;
4479  } else if (!value&&sparseThreshold_) {
4480    // delete copy
4481    sparseThreshold_=0;
4482    elementByRowL_.conditionalDelete();
4483    startRowL_.conditionalDelete();
4484    indexColumnL_.conditionalDelete();
4485    sparse_.conditionalDelete();
4486  } else if (value>0&&!sparseThreshold_) {
4487    if (value>1)
4488      sparseThreshold_=value;
4489    else
4490      sparseThreshold_=0;
4491    goSparse();
4492  }
4493#endif
4494}
4495void CoinAbcTypeFactorization::maximumPivots (  CoinSimplexInt value )
4496{
4497  if (value>0) {
4498    if (numberRows_)
4499      maximumMaximumPivots_=CoinMax(maximumMaximumPivots_,value);
4500    else
4501      maximumMaximumPivots_=value;
4502    maximumPivots_=value;
4503  }
4504  //printf("max %d maxmax %d\n",maximumPivots_,maximumMaximumPivots_);
4505}
4506void CoinAbcTypeFactorization::messageLevel (  CoinSimplexInt value )
4507{
4508  if (value>0&&value<16) {
4509    messageLevel_=value;
4510  }
4511}
4512// Reset all sparsity etc statistics
4513void CoinAbcTypeFactorization::resetStatistics()
4514{
4515#if ABC_SMALL<2
4516  setStatistics(true);
4517
4518  /// Below are all to collect
4519  ftranCountInput_=0.0;
4520  ftranCountAfterL_=0.0;
4521  ftranCountAfterR_=0.0;
4522  ftranCountAfterU_=0.0;
4523  ftranFTCountInput_=0.0;
4524  ftranFTCountAfterL_=0.0;
4525  ftranFTCountAfterR_=0.0;
4526  ftranFTCountAfterU_=0.0;
4527  btranCountInput_=0.0;
4528  btranCountAfterU_=0.0;
4529  btranCountAfterR_=0.0;
4530  btranCountAfterL_=0.0;
4531  ftranFullCountInput_=0.0;
4532  ftranFullCountAfterL_=0.0;
4533  ftranFullCountAfterR_=0.0;
4534  ftranFullCountAfterU_=0.0;
4535  btranFullCountInput_=0.0;
4536  btranFullCountAfterL_=0.0;
4537  btranFullCountAfterR_=0.0;
4538  btranFullCountAfterU_=0.0;
4539
4540  /// We can roll over factorizations
4541  numberFtranCounts_=0;
4542  numberFtranFTCounts_=0;
4543  numberBtranCounts_=0;
4544  numberFtranFullCounts_=0;
4545  numberFtranFullCounts_=0;
4546
4547  /// While these are averages collected over last
4548  ftranAverageAfterL_=INITIAL_AVERAGE;
4549  ftranAverageAfterR_=INITIAL_AVERAGE;
4550  ftranAverageAfterU_=INITIAL_AVERAGE;
4551  ftranFTAverageAfterL_=INITIAL_AVERAGE;
4552  ftranFTAverageAfterR_=INITIAL_AVERAGE;
4553  ftranFTAverageAfterU_=INITIAL_AVERAGE;
4554  btranAverageAfterU_=INITIAL_AVERAGE;
4555  btranAverageAfterR_=INITIAL_AVERAGE;
4556  btranAverageAfterL_=INITIAL_AVERAGE;
4557  ftranFullAverageAfterL_=INITIAL_AVERAGE;
4558  ftranFullAverageAfterR_=INITIAL_AVERAGE;
4559  ftranFullAverageAfterU_=INITIAL_AVERAGE;
4560  btranFullAverageAfterL_=INITIAL_AVERAGE;
4561  btranFullAverageAfterR_=INITIAL_AVERAGE;
4562  btranFullAverageAfterU_=INITIAL_AVERAGE;
4563#endif
4564}
4565// See if worth going sparse
4566void 
4567CoinAbcTypeFactorization::checkSparse()
4568{
4569#if ABC_SMALL<2
4570  // See if worth going sparse and when
4571  if (numberFtranCounts_>50) {
4572    ftranCountInput_= CoinMax(ftranCountInput_,1.0);
4573    ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,INITIAL_AVERAGE2);
4574    ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,INITIAL_AVERAGE2);
4575    ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,INITIAL_AVERAGE2);
4576    ftranFTCountInput_= CoinMax(ftranFTCountInput_,INITIAL_AVERAGE2);
4577    ftranFTAverageAfterL_ = CoinMax(ftranFTCountAfterL_/ftranFTCountInput_,INITIAL_AVERAGE2);
4578    ftranFTAverageAfterR_ = CoinMax(ftranFTCountAfterR_/ftranFTCountAfterL_,INITIAL_AVERAGE2);
4579    ftranFTAverageAfterU_ = CoinMax(ftranFTCountAfterU_/ftranFTCountAfterR_,INITIAL_AVERAGE2);
4580    if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
4581      btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,INITIAL_AVERAGE2);
4582      btranAverageAfterR_ = CoinMax(btranCountAfterR_/btranCountAfterU_,INITIAL_AVERAGE2);
4583      btranAverageAfterL_ = CoinMax(btranCountAfterL_/btranCountAfterR_,INITIAL_AVERAGE2);
4584    } else {
4585      // we have not done any useful btrans (values pass?)
4586      btranAverageAfterU_ = INITIAL_AVERAGE2;
4587      btranAverageAfterR_ = INITIAL_AVERAGE2;
4588      btranAverageAfterL_ = INITIAL_AVERAGE2;
4589    }
4590    ftranFullCountInput_= CoinMax(ftranFullCountInput_,1.0);
4591    ftranFullAverageAfterL_ = CoinMax(ftranFullCountAfterL_/ftranFullCountInput_,INITIAL_AVERAGE2);
4592    ftranFullAverageAfterR_ = CoinMax(ftranFullCountAfterR_/ftranFullCountAfterL_,INITIAL_AVERAGE2);
4593    ftranFullAverageAfterU_ = CoinMax(ftranFullCountAfterU_/ftranFullCountAfterR_,INITIAL_AVERAGE2);
4594    btranFullCountInput_= CoinMax(btranFullCountInput_,1.0);
4595    btranFullAverageAfterL_ = CoinMax(btranFullCountAfterL_/btranFullCountInput_,INITIAL_AVERAGE2);
4596    btranFullAverageAfterR_ = CoinMax(btranFullCountAfterR_/btranFullCountAfterL_,INITIAL_AVERAGE2);
4597    btranFullAverageAfterU_ = CoinMax(btranFullCountAfterU_/btranFullCountAfterR_,INITIAL_AVERAGE2);
4598  }
4599  // scale back
4600 
4601  ftranCountInput_ *= AVERAGE_SCALE_BACK;
4602  ftranCountAfterL_ *= AVERAGE_SCALE_BACK;
4603  ftranCountAfterR_ *= AVERAGE_SCALE_BACK;
4604  ftranCountAfterU_ *= AVERAGE_SCALE_BACK;
4605  ftranFTCountInput_ *= AVERAGE_SCALE_BACK;
4606  ftranFTCountAfterL_ *= AVERAGE_SCALE_BACK;
4607  ftranFTCountAfterR_ *= AVERAGE_SCALE_BACK;
4608  ftranFTCountAfterU_ *= AVERAGE_SCALE_BACK;
4609  btranCountInput_ *= AVERAGE_SCALE_BACK;
4610  btranCountAfterU_ *= AVERAGE_SCALE_BACK;
4611  btranCountAfterR_ *= AVERAGE_SCALE_BACK;
4612  btranCountAfterL_ *= AVERAGE_SCALE_BACK;
4613  ftranFullCountInput_ *= AVERAGE_SCALE_BACK;
4614  ftranFullCountAfterL_ *= AVERAGE_SCALE_BACK;
4615  ftranFullCountAfterR_ *= AVERAGE_SCALE_BACK;
4616  ftranFullCountAfterU_ *= AVERAGE_SCALE_BACK;
4617  btranFullCountInput_ *= AVERAGE_SCALE_BACK;
4618  btranFullCountAfterL_ *= AVERAGE_SCALE_BACK;
4619  btranFullCountAfterR_ *= AVERAGE_SCALE_BACK;
4620  btranFullCountAfterU_ *= AVERAGE_SCALE_BACK;
4621#endif
4622}
4623// Condition number - product of pivots after factorization
4624CoinSimplexDouble
4625CoinAbcTypeFactorization::conditionNumber() const
4626{
4627  CoinSimplexDouble condition = 1.0;
4628  const CoinFactorizationDouble * pivotRegion = pivotRegionAddress_;
4629  for (CoinSimplexInt i=0;i<numberRows_;i++) {
4630    condition *= pivotRegion[i];
4631  }
4632  condition = CoinMax(fabs(condition),1.0e-50);
4633  return 1.0/condition;
4634}
4635#ifndef NDEBUG
4636void 
4637CoinAbcTypeFactorization::checkMarkArrays() const
4638{
4639  CoinCheckZero * COIN_RESTRICT mark = markRowAddress_;
4640  for (CoinSimplexInt i=0;i<numberRows_;i++) {
4641    assert (!mark[i]);
4642  }
4643#if ABC_PARALLEL
4644  mark=reinterpret_cast<CoinCheckZero *>(reinterpret_cast<char *>(mark)+sizeSparseArray_);
4645  for (CoinSimplexInt i=0;i<numberRows_;i++) {
4646    assert (!mark[i]);
4647  }
4648  for (int i=1;i<FACTOR_CPU;i++) {
4649    for (CoinSimplexInt i=0;i<numberRows_;i++) {
4650      assert (!mark[i]);
4651    }
4652    mark=reinterpret_cast<CoinCheckZero *>(reinterpret_cast<char *>(mark)+sizeSparseArray_);
4653  }
4654#endif
4655}
4656#endif
4657#endif
4658#endif
Note: See TracBrowser for help on using the repository browser.