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

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

formatting

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