source: stable/1.15/Clp/src/CoinAbcBaseFactorization4.cpp @ 1995

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