source: trunk/Clp/src/CoinAbcBaseFactorization1.cpp @ 1878

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

minor changes to implement Aboca

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