source: stable/1.15/Clp/src/CoinAbcBaseFactorization1.cpp @ 2006

Last change on this file since 2006 was 2006, checked in by forrest, 6 years ago

changes for parallel and idiot

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