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

Last change on this file since 1910 was 1910, checked in by stefan, 7 years ago
  • add configure option --enable-aboca={1,2,3,4,yes,no}
  • compile Aboca source only if --enable-aboca set (instead of compiling empty source files)
  • fix svn properties
  • Property svn:keywords set to Id
File size: 154.5 KB
Line 
1/* $Id: CoinAbcBaseFactorization1.cpp 1910 2013-01-27 02:00:13Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5#ifdef ABC_JUST_ONE_FACTORIZATION
6#include "CoinAbcCommonFactorization.hpp"
7#define CoinAbcTypeFactorization CoinAbcBaseFactorization
8#define ABC_SMALL -1
9#include "CoinAbcBaseFactorization.hpp"
10#endif
11#ifdef CoinAbcTypeFactorization
12#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    for (int iTry=numberZeroOther-1;iTry>=0;iTry--) {
1486      pivotWhile(0,numberInPivotRow,numberInPivotColumn,lengthArea,giveUp,
1487                 eArea,multipliersL);
1488      // find new pivot row
1489      int jPivotColumn=others[iTry];
1490      CoinFactorizationDouble * COIN_RESTRICT area =eArea+jPivotColumn*lengthArea;
1491      double largest=tolerance;
1492      int jPivotRow=-1;
1493      for (CoinSimplexInt j = 0; j < numberInPivotColumn; j++ ) {
1494        CoinFactorizationDouble value = area[j];
1495        CoinSimplexDouble absValue = fabs ( value );
1496        if (absValue>largest) {
1497          largest=absValue;
1498          jPivotRow=j;
1499        }
1500      }
1501      if (jPivotRow>=0) {
1502        markRow[pivotRow] = LARGE_UNSET;
1503        //modify linked list for pivots
1504        deleteLink ( pivotRow );
1505        deleteLink ( pivotColumn + numberRows_ );
1506        // put away last
1507        afterPivot(pivotRow,pivotColumn);
1508        // new ones
1509        pivotRow=indexL[jPivotRow];
1510        pivotColumn=saveColumn[jPivotColumn];
1511        // now do stuff for pivot
1512        realPivotRow=fromSmallToBigRow_[pivotRow];
1513        //int realPivotColumn=fromSmallToBigColumn[pivotColumn];
1514        //store pivot columns (so can easily compress)
1515        CoinBigIndex startColumn = startColumnU[pivotColumn];
1516        CoinBigIndex endColumn = startColumn + tempColumnCount[jPivotColumn];
1517        CoinSimplexInt next = nextRow[pivotRow];
1518        CoinSimplexInt last = lastRow[pivotRow];
1519
1520        nextRow[last] = next;
1521        lastRow[next] = last;
1522#ifdef SMALL_PERMUTE
1523        permuteAddress_[realPivotRow] = numberGoodU_;   //use for permute
1524#else
1525        permuteAddress_[pivotRow] = numberGoodU_;       //use for permute
1526#endif
1527        lastRow[pivotRow] = -2;
1528        numberInRow[pivotRow] = 0;
1529        //store column in L, compress in U and take column out
1530        CoinBigIndex l = lengthL_;
1531        if ( l + numberInPivotColumn > lengthAreaL_ ) {
1532          //need more memory
1533          if ((messageLevel_&4)!=0) 
1534            printf("more memory needed in middle of invert\n");
1535          return -99;
1536        }
1537        //CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnL_.array();
1538        numberGoodL_++;
1539        // move to L and move last
1540        numberInPivotColumn--;
1541        if (jPivotRow!=numberInPivotColumn) {
1542          // swap
1543          for (int jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1544            CoinFactorizationDouble * COIN_RESTRICT area =eArea+jColumn*lengthArea;
1545            CoinFactorizationDouble thisValue = area[jPivotRow];
1546            area[jPivotRow]=area[numberInPivotColumn];
1547            area[numberInPivotColumn]=thisValue;
1548          }
1549          // do indexL in last one
1550          int tempI=indexL[numberInPivotColumn];
1551          CoinFactorizationDouble tempD=multipliersL[numberInPivotColumn];
1552          indexL[numberInPivotColumn]=indexL[jPivotRow];
1553          multipliersL[numberInPivotColumn]=multipliersL[jPivotRow];
1554          indexL[jPivotRow]=tempI;
1555          multipliersL[jPivotRow]=tempD;
1556          jPivotRow=numberInPivotColumn;
1557        }
1558        numberInPivotRow--;
1559        CoinFactorizationDouble pivotElement = area[numberInPivotColumn];
1560        area[numberInPivotColumn]=0.0;
1561        // put last one in now
1562        CoinFactorizationDouble * COIN_RESTRICT areaLast =eArea+numberInPivotRow*lengthArea;
1563        // swap columns
1564        saveColumn[jPivotColumn]=saveColumn[numberInPivotRow];
1565        saveColumn[numberInPivotRow]=pivotColumn;
1566        //int temp=tempColumnCount[jPivotColumn];
1567        tempColumnCount[jPivotColumn]=tempColumnCount[numberInPivotRow];
1568        totalElements_-=numberInColumn[pivotColumn];
1569        mark[pivotColumn]=0;
1570        for (int jRow=0;jRow<numberInPivotColumn+1;jRow++) {
1571          CoinFactorizationDouble temp=areaLast[jRow];
1572          areaLast[jRow]=area[jRow];
1573          area[jRow]=temp;
1574        }
1575        areaLast[numberInPivotColumn]=0.0;
1576        // swap rows in area and save
1577        for (int jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1578          CoinFactorizationDouble * COIN_RESTRICT area =eArea+jColumn*lengthArea;
1579          CoinFactorizationDouble thisPivotValue = area[numberInPivotColumn];
1580          area[numberInPivotColumn]=thisPivotValue;
1581          if (fabs(thisPivotValue)>tolerance) {
1582            CoinSimplexInt iColumn = saveColumn[jColumn];
1583            CoinBigIndex startColumn = startColumnU[iColumn];
1584            int nKeep=tempColumnCount[jColumn];
1585            if (nKeep) {
1586              if (nKeep>1) {
1587                // need to move one to end
1588                elementU[startColumn+nKeep]=elementU[startColumn+1];
1589                indexRowU[startColumn+nKeep]=indexRowU[startColumn+1];
1590              }
1591              // move largest
1592              elementU[startColumn+1]=elementU[startColumn];
1593              indexRowU[startColumn+1]=indexRowU[startColumn];
1594            }
1595            elementU[startColumn] = thisPivotValue;
1596            indexRowU[startColumn] = realPivotRow;
1597            startColumn++;
1598            startColumnU[iColumn]=startColumn;
1599            numberInColumnPlus[iColumn]++;
1600            //numberInColumn[iColumn]--;
1601          }
1602        }
1603        // put away last (but remember count is one out)
1604        CoinAbcMemcpy(elementL+lSave,multipliersL,numberInPivotColumn+1);
1605        lSave=l;
1606        for (int i=0;i<numberInPivotColumn;i++) {
1607          CoinFactorizationDouble value=areaLast[i];
1608          areaLast[i]=0.0;
1609          int iRow=indexL[i];
1610          indexRowL[l] = iRow;
1611          assert (iRow>=0);
1612          multipliersL[l-lSave] = value;
1613          l++;
1614        }
1615        startColumnL[numberGoodL_] = l;
1616        lengthL_ = l;
1617        CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1618        pivotRegionAddress_[numberGoodU_] = pivotMultiplier;
1619        numberInColumn[pivotColumn] = 0;
1620        indexL = &indexRowL[lSave];
1621        //adjust
1622        for (CoinSimplexInt j = 0; j < numberInPivotColumn; j++ ) {
1623          multipliersL[j] *= pivotMultiplier;
1624        }
1625        // zero out last
1626        multipliersL[numberInPivotColumn]=0.0;
1627      } else {
1628        // singular
1629        // give up
1630        break;
1631      }
1632    }
1633#if 0
1634    // zero out extra bits
1635    int add1=numberInPivotColumn>>5;
1636    unsigned int * aBitsThis = aBits;
1637    int left=numberInPivotColumn-(add1<<5);
1638    unsigned int mark1 = (1<<(left-1)-1);
1639
1640    int jColumn;
1641    for (jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1642      memset(aBitsThis+add1,0,(intsPerColumn-add1)*sizeof(int));
1643      aBitsThis += intsPerColumn;
1644    }
1645    for (; jColumn < numberInPivotRow+numberZeroOther; jColumn++ ) {
1646      memset(aBitsThis,0,intsPerColumn*sizeof(int));
1647      aBitsThis += intsPerColumn;
1648    }
1649#endif
1650    // last bit
1651    inAreaAlready=true;
1652  }
1653#endif
1654#if 0
1655  for (int i=0;i<numberRowsSmall_;i++)
1656    assert (numberInColumn[i]>=0);
1657#endif
1658  //if (added<1000)
1659  //giveUp=CoinMax(); //??
1660#ifdef DENSE_TRY
1661  if (!inAreaAlready) {
1662#endif
1663    pivotSome(0,numberInPivotRow,numberInPivotColumn,lengthArea,giveUp,
1664              eArea,saveColumn,startColumnU,tempColumnCount,
1665              elementU,numberInColumn,indexRowU,aBits,indexL,
1666              multipliersL,tolerance);
1667#ifdef DENSE_TRY
1668  } else {
1669    pivotSomeAfter(0,numberInPivotRow,numberInPivotColumn,lengthArea,giveUp,
1670              eArea,saveColumn,startColumnU,tempColumnCount,
1671              elementU,numberInColumn,indexRowU,aBits,indexL,
1672              multipliersL,tolerance);
1673  }
1674#endif
1675#if 0
1676  for (int i=0;i<numberRowsSmall_;i++)
1677    assert (numberInColumn[i]>=0);
1678#endif
1679  //serial
1680  for (int jRow = 0; jRow < numberInPivotColumn; jRow++ ) {
1681    CoinSimplexInt iRow = indexL[jRow];
1682    CoinSimplexInt next = nextRow[iRow];
1683    CoinBigIndex space = startRowU[next] - startRowU[iRow];
1684    // assume full
1685    int numberNeeded = numberInPivotRow+numberInRow[iRow];
1686    tempCount[jRow]=numberInRow[iRow];
1687    numberInRow[iRow]=numberNeeded;
1688    if ( space < numberNeeded ) {
1689      if ( !getRowSpace ( iRow, numberNeeded ) ) {
1690        return -99;
1691      }
1692    }
1693    lastEntryByRowU_ = CoinMax(startRowU[iRow]+numberNeeded,lastEntryByRowU_);
1694    // in case I got it wrong!
1695    if (lastEntryByRowU_>lengthAreaU_) {
1696      return -99;
1697    }
1698  }
1699  // move counts back
1700  for (int jRow = 0; jRow < numberInPivotColumn; jRow++ ) {
1701    CoinSimplexInt iRow = indexL[jRow];
1702    numberInRow[iRow]=tempCount[jRow];
1703  }
1704  // parallel
1705  cilk_for_query (int jRow = 0; jRow < numberInPivotColumn; jRow++ ) {
1706    CoinSimplexInt iRow = indexL[jRow];
1707    CoinBigIndex startRow = startRowU[iRow];
1708    CoinBigIndex put = startRow + numberInRow[iRow];
1709    int iWord = jRow>>5;
1710    unsigned int setBit=1<<(jRow&31);
1711    unsigned int * aBitsThis = aBits+iWord;
1712    for (int i=0;i<numberInPivotRow;i++) {
1713      if ((aBitsThis[0]&setBit)!=0) {
1714        indexColumnU[put++]=saveColumn[i];
1715      }
1716      aBitsThis += intsPerColumn;
1717    }
1718    markRow[iRow] = LARGE_UNSET;
1719    numberInRow[iRow] = put-startRow;
1720  }
1721  //memset(aBits,0x00,intsPerColumn*4*(numberInPivotRow+numberZeroOther));
1722  added=0;
1723  //int addedX=0;
1724  for (int jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1725    CoinSimplexInt iColumn = saveColumn[jColumn];
1726    for (int i=startColumnU[iColumn];i<startColumnU[iColumn]+numberInColumn[iColumn];i++)
1727      assert (indexRowU[i]>=0&&indexRowU[i]<numberRowsSmall_);
1728    mark[iColumn]=0;
1729    added += tempColumnCount[jColumn];
1730    modifyLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1731    //addedX += numberInColumn[iColumn]-tempxx[jColumn];
1732  }
1733#ifndef NDEBUG
1734  for (int i=0;i<added;i++)
1735    assert (!eArea[i]);
1736  for (int i=0;i<firstZeroed_;i++)
1737    assert (!elementUColumnPlusAddress_[lengthAreaUPlus_-firstZeroed_+i]);
1738  for (int i=0;i<numberRowsSmall_;i++)
1739    assert (!mark[i]);
1740#endif
1741  //double ratio=((double)addedX)/((double)(numberInPivotRow*numberInPivotColumn));
1742#if 0 //def DENSE_TRY
1743  if (numberZeroOther>10000)
1744    printf("INrow %d incolumn %d nzero %d\n",
1745           numberInPivotRow,numberInPivotColumn,numberZeroOther);
1746#endif
1747  //printf("INrow %d incolumn %d nzero %d ratio %g\n",
1748  //     numberInPivotRow,numberInPivotColumn,numberZeroOther,ratio);
1749  for (int i=0;i<numberInPivotColumn;i++) {
1750    int iRow=indexL[i];
1751    modifyLink(iRow,numberInRow[iRow]);
1752  }
1753  CoinAbcMemcpy(elementL+lSave,multipliersL,numberInPivotColumn);
1754#ifdef DENSE_TRY
1755  int iLBad=saveGoodL;
1756  put=startColumnL[iLBad];
1757  for (int iL=iLBad;iL<numberGoodL_;iL++) {
1758    CoinBigIndex start=startColumnL[iL];
1759    startColumnL[iL]=put;
1760    for (CoinBigIndex j=start;j<startColumnL[iL+1];j++) {
1761      if (elementL[j]) {
1762        elementL[put]=elementL[j];
1763        indexRowL[put++]=fromSmallToBigRow_[indexRowL[j]];
1764      }
1765    }
1766  }
1767  startColumnL[numberGoodL_]=put;
1768  lengthL_=put;
1769#else
1770  // redo L indices
1771  for (CoinBigIndex j=lSave;j<startColumnL[numberGoodL_];j++) {
1772    //assert (fabs(elementL[j])>=tolerance) ;
1773    int iRow=indexRowL[j];
1774    indexRowL[j] = fromSmallToBigRow_[iRow];
1775  }
1776#endif
1777  markRow[pivotRow] = LARGE_UNSET;
1778  //modify linked list for pivots
1779  deleteLink ( pivotRow );
1780  deleteLink ( pivotColumn + numberRows_ );
1781  totalElements_ += added;
1782  return status;
1783}
1784#endif
1785// nonparallel pivot
1786int
1787CoinAbcTypeFactorization::pivot ( CoinSimplexInt pivotRow,
1788                                  CoinSimplexInt pivotColumn,
1789                                  CoinBigIndex pivotRowPosition,
1790                                  CoinBigIndex pivotColumnPosition,
1791                                  CoinFactorizationDouble * COIN_RESTRICT work,
1792                                  CoinSimplexUnsignedInt * COIN_RESTRICT workArea2,
1793                                  CoinSimplexInt increment2,
1794                                  int * COIN_RESTRICT markRow)
1795{
1796  CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnU_.array();
1797  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnU_.array();
1798  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumn_.array();
1799  CoinFactorizationDouble * COIN_RESTRICT elementU = elementU_.array();
1800  CoinSimplexInt * COIN_RESTRICT indexRowU = indexRowU_.array();
1801  CoinBigIndex * COIN_RESTRICT startRowU = startRowU_.array();
1802  CoinSimplexInt * COIN_RESTRICT numberInRow = numberInRow_.array();
1803  CoinFactorizationDouble * COIN_RESTRICT elementL = elementL_.array();
1804  CoinSimplexInt * COIN_RESTRICT indexRowL = indexRowL_.array();
1805  CoinSimplexInt * COIN_RESTRICT saveColumn = saveColumn_.array();
1806  CoinSimplexInt * COIN_RESTRICT nextRow = nextRow_.array();
1807  CoinSimplexInt * COIN_RESTRICT lastRow = lastRow_.array() ;
1808#ifdef SMALL_PERMUTE
1809  int realPivotRow=fromSmallToBigRow_[pivotRow];
1810  //int realPivotColumn=fromSmallToBigColumn[pivotColumn];
1811#endif
1812  //store pivot columns (so can easily compress)
1813  CoinSimplexInt numberInPivotRow = numberInRow[pivotRow] - 1;
1814  CoinBigIndex startColumn = startColumnU[pivotColumn];
1815  CoinSimplexInt numberInPivotColumn = numberInColumn[pivotColumn] - 1;
1816  // temp
1817  CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
1818  CoinBigIndex put = 0;
1819  CoinBigIndex startRow = startRowU[pivotRow];
1820  CoinBigIndex endRow = startRow + numberInPivotRow + 1;
1821#if ABC_SMALL<2
1822  // temp - switch off marker for valid row/column lookup
1823  sparseThreshold_=-1;
1824#endif
1825
1826  if ( pivotColumnPosition < 0 ) {
1827    for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1828      CoinSimplexInt iColumn = indexColumnU[pivotColumnPosition];
1829      if ( iColumn != pivotColumn ) {
1830        saveColumn[put++] = iColumn;
1831      } else {
1832#if BOTH_WAYS
1833#endif
1834        break;
1835      }
1836    }
1837  } else {
1838    for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
1839      saveColumn[put++] = indexColumnU[i];
1840    }
1841  }
1842  assert (pivotColumnPosition<endRow);
1843  assert (indexColumnU[pivotColumnPosition]==pivotColumn);
1844  pivotColumnPosition++;
1845  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1846    saveColumn[put++] = indexColumnU[pivotColumnPosition];
1847  }
1848  //take out this bit of indexColumnU
1849  CoinSimplexInt next = nextRow[pivotRow];
1850  CoinSimplexInt last = lastRow[pivotRow];
1851
1852  nextRow[last] = next;
1853  lastRow[next] = last;
1854#ifdef SMALL_PERMUTE
1855  permuteAddress_[realPivotRow] = numberGoodU_; //use for permute
1856#else
1857  permuteAddress_[pivotRow] = numberGoodU_;     //use for permute
1858#endif
1859  lastRow[pivotRow] = -2;
1860  numberInRow[pivotRow] = 0;
1861  //store column in L, compress in U and take column out
1862  CoinBigIndex l = lengthL_;
1863
1864  if ( l + numberInPivotColumn > lengthAreaL_ ) {
1865    //need more memory
1866    if ((messageLevel_&4)!=0) 
1867      printf("more memory needed in middle of invert\n");
1868    return -99;
1869  }
1870  //l+=currentAreaL_->elementByColumn-elementL;
1871  CoinBigIndex lSave = l;
1872
1873  CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnL_.array();
1874  startColumnL[numberGoodL_] = l;       //for luck and first time
1875  numberGoodL_++;
1876  startColumnL[numberGoodL_] = l + numberInPivotColumn;
1877  lengthL_ += numberInPivotColumn;
1878  if ( pivotRowPosition < 0 ) {
1879    for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1880      CoinSimplexInt iRow = indexRowU[pivotRowPosition];
1881      if ( iRow != pivotRow ) {
1882        indexRowL[l] = iRow;
1883        elementL[l] = elementU[pivotRowPosition];
1884        markRow[iRow] = l - lSave;
1885        l++;
1886        //take out of row list
1887        CoinBigIndex start = startRowU[iRow];
1888        CoinBigIndex end = start + numberInRow[iRow];
1889        CoinBigIndex where = start;
1890
1891        while ( indexColumnU[where] != pivotColumn ) {
1892          where++;
1893        }                       /* endwhile */
1894#if DEBUG_COIN
1895        if ( where >= end ) {
1896          abort (  );
1897        }
1898#endif
1899#if BOTH_WAYS
1900#endif
1901        indexColumnU[where] = indexColumnU[end - 1];
1902        numberInRow[iRow]--;
1903      } else {
1904        break;
1905      }
1906    }
1907  } else {
1908    CoinBigIndex i;
1909
1910    for ( i = startColumn; i < pivotRowPosition; i++ ) {
1911      CoinSimplexInt iRow = indexRowU[i];
1912
1913      markRow[iRow] = l - lSave;
1914      indexRowL[l] = iRow;
1915      elementL[l] = elementU[i];
1916      l++;
1917      //take out of row list
1918      CoinBigIndex start = startRowU[iRow];
1919      CoinBigIndex end = start + numberInRow[iRow];
1920      CoinBigIndex where = start;
1921
1922      while ( indexColumnU[where] != pivotColumn ) {
1923        where++;
1924      }                         /* endwhile */
1925#if DEBUG_COIN
1926      if ( where >= end ) {
1927        abort (  );
1928      }
1929#endif
1930#if BOTH_WAYS
1931#endif
1932      indexColumnU[where] = indexColumnU[end - 1];
1933      numberInRow[iRow]--;
1934      assert (numberInRow[iRow]>=0);
1935    }
1936  }
1937  assert (pivotRowPosition<endColumn);
1938  assert (indexRowU[pivotRowPosition]==pivotRow);
1939  CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
1940  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1941
1942  pivotRegionAddress_[numberGoodU_] = pivotMultiplier;
1943  pivotRowPosition++;
1944  for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1945    CoinSimplexInt iRow = indexRowU[pivotRowPosition];
1946   
1947    markRow[iRow] = l - lSave;
1948    indexRowL[l] = iRow;
1949    elementL[l] = elementU[pivotRowPosition];
1950    l++;
1951    //take out of row list
1952    CoinBigIndex start = startRowU[iRow];
1953    CoinBigIndex end = start + numberInRow[iRow];
1954    CoinBigIndex where = start;
1955    // could vectorize?
1956    while ( indexColumnU[where] != pivotColumn ) {
1957      where++;
1958    }                           /* endwhile */
1959#if DEBUG_COIN
1960    if ( where >= end ) {
1961      abort (  );
1962    }
1963#endif
1964#if BOTH_WAYS
1965#endif
1966    indexColumnU[where] = indexColumnU[end - 1];
1967    numberInRow[iRow]--;
1968    assert (numberInRow[iRow]>=0);
1969  }
1970  markRow[pivotRow] = LARGE_SET;
1971  //compress pivot column (move pivot to front including saved)
1972  numberInColumn[pivotColumn] = 0;
1973  //use end of L for temporary space
1974  CoinSimplexInt * COIN_RESTRICT indexL = &indexRowL[lSave];
1975  CoinFactorizationDouble * COIN_RESTRICT multipliersL = &elementL[lSave];
1976
1977  //adjust
1978  CoinSimplexInt j;
1979
1980  for ( j = 0; j < numberInPivotColumn; j++ ) {
1981    multipliersL[j] *= pivotMultiplier;
1982  }
1983  CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
1984  //zero out fill
1985  CoinBigIndex iErase;
1986  for ( iErase = 0; iErase < increment2 * numberInPivotRow;
1987        iErase++ ) {
1988    workArea2[iErase] = 0;
1989  }
1990  CoinSimplexUnsignedInt * COIN_RESTRICT temp2 = workArea2;
1991  CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumn_.array();
1992
1993  //pack down and move to work
1994  CoinSimplexInt jColumn;
1995  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1996    CoinSimplexInt iColumn = saveColumn[jColumn];
1997    CoinBigIndex startColumn = startColumnU[iColumn];
1998    CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
1999    CoinSimplexInt iRow = indexRowU[startColumn];
2000    CoinFactorizationDouble value = elementU[startColumn];
2001    CoinSimplexDouble largest;
2002    CoinBigIndex put = startColumn;
2003    CoinBigIndex positionLargest = -1;
2004    CoinFactorizationDouble thisPivotValue = 0.0;
2005
2006    //compress column and find largest not updated
2007    bool checkLargest;
2008    CoinSimplexInt mark = markRow[iRow];
2009
2010    if ( mark == LARGE_UNSET ) {
2011      largest = fabs ( value );
2012      positionLargest = put;
2013      put++;
2014      checkLargest = false;
2015    } else {
2016      //need to find largest
2017      largest = 0.0;
2018      checkLargest = true;
2019      if ( mark != LARGE_SET ) {
2020        //will be updated
2021        work[mark] = value;
2022        CoinSimplexInt word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
2023        CoinSimplexInt bit = mark & COINFACTORIZATION_MASK_PER_INT;
2024
2025        temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
2026        added--;
2027      } else {
2028        thisPivotValue = value;
2029      }
2030    }
2031    CoinBigIndex i;
2032    for ( i = startColumn + 1; i < endColumn; i++ ) {
2033      iRow = indexRowU[i];
2034      value = elementU[i];
2035      CoinSimplexInt mark = markRow[iRow];
2036
2037      if ( mark == LARGE_UNSET ) {
2038        //keep
2039        indexRowU[put] = iRow;
2040#if BOTH_WAYS
2041#endif
2042        elementU[put] = value;
2043        if ( checkLargest ) {
2044          CoinSimplexDouble absValue = fabs ( value );
2045
2046          if ( absValue > largest ) {
2047            largest = absValue;
2048            positionLargest = put;
2049          }
2050        }
2051        put++;
2052      } else if ( mark != LARGE_SET ) {
2053        //will be updated
2054        work[mark] = value;
2055        CoinSimplexInt word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
2056        CoinSimplexInt bit = mark & COINFACTORIZATION_MASK_PER_INT;
2057
2058        temp2[word] = temp2[word] | ( 1 << bit );       //say already in counts
2059        added--;
2060      } else {
2061        thisPivotValue = value;
2062      }
2063    }
2064    //slot in pivot
2065    elementU[put] = elementU[startColumn];
2066    indexRowU[put] = indexRowU[startColumn];
2067#if BOTH_WAYS
2068#endif
2069    if ( positionLargest == startColumn ) {
2070      positionLargest = put;    //follow if was largest
2071    }
2072    put++;
2073    elementU[startColumn] = thisPivotValue;
2074#ifdef SMALL_PERMUTE
2075    indexRowU[startColumn] = realPivotRow;
2076#else
2077    indexRowU[startColumn] = pivotRow;
2078#endif
2079#if BOTH_WAYS
2080#endif
2081    //clean up counts
2082    startColumn++;
2083    numberInColumn[iColumn] = put - startColumn;
2084    CoinSimplexInt *  COIN_RESTRICT numberInColumnPlus = numberInColumnPlus_.array();
2085    numberInColumnPlus[iColumn]++;
2086    startColumnU[iColumn]++;
2087    //how much space have we got
2088    CoinSimplexInt next = nextColumn[iColumn];
2089    CoinBigIndex space;
2090
2091    space = startColumnU[next] - put - numberInColumnPlus[next];
2092    //assume no zero elements
2093    if ( numberInPivotColumn > space ) {
2094      //getColumnSpace also moves fixed part
2095      if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
2096        return -99;
2097      }
2098      //redo starts
2099      positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
2100      startColumn = startColumnU[iColumn];
2101      put = startColumn + numberInColumn[iColumn];
2102    }
2103    CoinSimplexDouble tolerance = zeroTolerance_;
2104
2105    //CoinSimplexInt * COIN_RESTRICT nextCount = this->nextCount();
2106    for ( j = 0; j < numberInPivotColumn; j++ ) {
2107      value = work[j] - thisPivotValue * multipliersL[j];
2108      CoinSimplexDouble absValue = fabs ( value );
2109
2110      if ( absValue > tolerance ) {
2111        work[j] = 0.0;
2112        assert (put<lengthAreaU_); 
2113        elementU[put] = value;
2114        indexRowU[put] = indexL[j];
2115#if BOTH_WAYS
2116#endif
2117        if ( absValue > largest ) {
2118          largest = absValue;
2119          positionLargest = put;
2120        }
2121        put++;
2122      } else {
2123        work[j] = 0.0;
2124        added--;
2125        CoinSimplexInt word = j >> COINFACTORIZATION_SHIFT_PER_INT;
2126        CoinSimplexInt bit = j & COINFACTORIZATION_MASK_PER_INT;
2127
2128        if ( temp2[word] & ( 1 << bit ) ) {
2129          //take out of row list
2130          iRow = indexL[j];
2131          CoinBigIndex start = startRowU[iRow];
2132          CoinBigIndex end = start + numberInRow[iRow];
2133          CoinBigIndex where = start;
2134
2135          while ( indexColumnU[where] != iColumn ) {
2136            where++;
2137          }                     /* endwhile */
2138#if DEBUG_COIN
2139          if ( where >= end ) {
2140            abort (  );
2141          }
2142#endif
2143#if BOTH_WAYS
2144#endif
2145          indexColumnU[where] = indexColumnU[end - 1];
2146          numberInRow[iRow]--;
2147        } else {
2148          //make sure won't be added
2149          CoinSimplexInt word = j >> COINFACTORIZATION_SHIFT_PER_INT;
2150          CoinSimplexInt bit = j & COINFACTORIZATION_MASK_PER_INT;
2151
2152          temp2[word] = temp2[word] | ( 1 << bit );     //say already in counts
2153        }
2154      }
2155    }
2156    numberInColumn[iColumn] = put - startColumn;
2157    //move largest
2158    if ( positionLargest >= 0 ) {
2159      value = elementU[positionLargest];
2160      iRow = indexRowU[positionLargest];
2161      elementU[positionLargest] = elementU[startColumn];
2162      indexRowU[positionLargest] = indexRowU[startColumn];
2163      elementU[startColumn] = value;
2164      indexRowU[startColumn] = iRow;
2165#if BOTH_WAYS
2166#endif
2167    }
2168    //linked list for column
2169    //if ( nextCount[iColumn + numberRows_] != -2 ) {
2170      //modify linked list
2171      //deleteLink ( iColumn + numberRows_ );
2172      //modifyLink ( iColumn + numberRows_, numberInColumn[iColumn] );
2173      //}
2174    temp2 += increment2;
2175  }
2176  //get space for row list
2177  CoinSimplexUnsignedInt * COIN_RESTRICT putBase = workArea2;
2178  CoinSimplexInt bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
2179  CoinSimplexInt i = 0;
2180
2181  // do linked lists and update counts
2182  while ( bigLoops ) {
2183    bigLoops--;
2184    CoinSimplexInt bit;
2185    for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
2186      CoinSimplexUnsignedInt *putThis = putBase;
2187      CoinSimplexInt iRow = indexL[i];
2188
2189      //get space
2190      CoinSimplexInt number = 0;
2191      CoinSimplexInt jColumn;
2192
2193      for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
2194        CoinSimplexUnsignedInt test = *putThis;
2195
2196        putThis += increment2;
2197        test = 1 - ( ( test >> bit ) & 1 );
2198        number += test;
2199      }
2200      CoinSimplexInt next = nextRow[iRow];
2201      CoinBigIndex space;
2202
2203      space = startRowU[next] - startRowU[iRow];
2204      number += numberInRow[iRow];
2205      if ( space < number ) {
2206        if ( !getRowSpace ( iRow, number ) ) {
2207          return -99;
2208        }
2209      }
2210      // now do
2211      putThis = putBase;
2212      next = nextRow[iRow];
2213      number = numberInRow[iRow];
2214      CoinBigIndex end = startRowU[iRow] + number;
2215      CoinSimplexInt saveIndex = indexColumnU[startRowU[next]];
2216
2217      //add in
2218      for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
2219        CoinSimplexUnsignedInt test = *putThis;
2220
2221        putThis += increment2;
2222        test = 1 - ( ( test >> bit ) & 1 );
2223        indexColumnU[end] = saveColumn[jColumn];
2224#if BOTH_WAYS
2225#endif
2226        end += test;
2227      }
2228      lastEntryByRowU_ = CoinMax(end,lastEntryByRowU_);
2229      assert (lastEntryByRowU_<=lengthAreaU_);
2230      //put back next one in case zapped
2231      indexColumnU[startRowU[next]] = saveIndex;
2232#if BOTH_WAYS
2233#endif
2234      markRow[iRow] = LARGE_UNSET;
2235      number = end - startRowU[iRow];
2236      numberInRow[iRow] = number;
2237      //deleteLink ( iRow );
2238      //modifyLink ( iRow, number );
2239    }
2240    putBase++;
2241  }                             /* endwhile */
2242  CoinSimplexInt bit;
2243
2244  for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
2245    CoinSimplexUnsignedInt *putThis = putBase;
2246    CoinSimplexInt iRow = indexL[i];
2247
2248    //get space
2249    CoinSimplexInt number = 0;
2250    CoinSimplexInt jColumn;
2251
2252    for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
2253      CoinSimplexUnsignedInt test = *putThis;
2254
2255      putThis += increment2;
2256      test = 1 - ( ( test >> bit ) & 1 );
2257      number += test;
2258    }
2259    CoinSimplexInt next = nextRow[iRow];
2260    CoinBigIndex space;
2261
2262    space = startRowU[next] - startRowU[iRow];
2263    number += numberInRow[iRow];
2264    if ( space < number ) {
2265      if ( !getRowSpace ( iRow, number ) ) {
2266        return -99;
2267      }
2268    }
2269    // now do
2270    putThis = putBase;
2271    next = nextRow[iRow];
2272    number = numberInRow[iRow];
2273    CoinBigIndex end = startRowU[iRow] + number;
2274    CoinSimplexInt saveIndex;
2275
2276    saveIndex = indexColumnU[startRowU[next]];
2277    //if (numberRowsSmall_==2791&&iRow==1160)
2278    //printf("start %d end %d\n",startRowU[iRow],end);
2279    //add in
2280    for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
2281      CoinSimplexUnsignedInt test = *putThis;
2282
2283      putThis += increment2;
2284      test = 1 - ( ( test >> bit ) & 1 );
2285
2286      indexColumnU[end] = saveColumn[jColumn];
2287#if BOTH_WAYS
2288#endif
2289      end += test;
2290      //if (numberRowsSmall_==2791&&iRow==1160)
2291      //printf("start %d end %d test %d col %d\n",startRowU[iRow],end,test,saveColumn[jColumn]);
2292    }
2293    indexColumnU[startRowU[next]] = saveIndex;
2294    lastEntryByRowU_ = CoinMax(end,lastEntryByRowU_);
2295    assert (lastEntryByRowU_<=lengthAreaU_);
2296#if BOTH_WAYS
2297#endif
2298    markRow[iRow] = LARGE_UNSET;
2299    number = end - startRowU[iRow];
2300    numberInRow[iRow] = number;
2301    //deleteLink ( iRow );
2302    //modifyLink ( iRow, number );
2303  }
2304  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
2305    CoinSimplexInt iColumn = saveColumn[jColumn];
2306    modifyLink ( iColumn + numberRows_, numberInColumn[iColumn] );
2307  }
2308  for (int i=0;i<numberInPivotColumn;i++) {
2309    int iRow=indexL[i];
2310    modifyLink(iRow,numberInRow[iRow]);
2311  }
2312#ifdef SMALL_PERMUTE
2313  // redo L indices
2314  for (CoinBigIndex j=lSave;j<startColumnL[numberGoodL_];j++) {
2315    int iRow=indexRowL[j];
2316    indexRowL[j] = fromSmallToBigRow_[iRow];
2317  }
2318#endif
2319  markRow[pivotRow] = LARGE_UNSET;
2320  //modify linked list for pivots
2321  deleteLink ( pivotRow );
2322  deleteLink ( pivotColumn + numberRows_ );
2323  totalElements_ += added;
2324  return 0;
2325}
2326// After pivoting
2327void
2328CoinAbcTypeFactorization::afterPivot( CoinSimplexInt pivotRow,
2329                                      CoinSimplexInt pivotColumn )
2330{
2331#ifdef SMALL_PERMUTE
2332  int realPivotRow=fromSmallToBigRow_[pivotRow];
2333  int realPivotColumn=fromSmallToBigColumn_[pivotColumn];
2334  assert (permuteAddress_[realPivotRow]==numberGoodU_);
2335  pivotColumnAddress_[numberGoodU_] = realPivotColumn;
2336#else
2337  assert (permuteAddress_[pivotRow]==numberGoodU_);
2338  pivotColumnAddress_[numberGoodU_] = pivotColumn;
2339#endif
2340#ifdef ABC_USE_FUNCTION_POINTERS
2341  int number = numberInColumnPlusAddress_[pivotColumn];
2342  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
2343#ifdef SMALL_PERMUTE
2344#if ABC_USE_FUNCTION_POINTERS
2345  if (number<9) {
2346    scatter[realPivotRow].functionPointer=AbcScatterLowSubtract[number];
2347  } else {
2348    scatter[realPivotRow].functionPointer=AbcScatterHighSubtract[number&3];
2349  }
2350#endif
2351  scatter[realPivotRow].offset=lastEntryByColumnUPlus_;
2352  scatter[realPivotRow].number=number;
2353#else
2354#if ABC_USE_FUNCTION_POINTERS
2355  if (number<9) {
2356    scatter[pivotRow].functionPointer=AbcScatterLowSubtract[number];
2357  } else {
2358    scatter[pivotRow].functionPointer=AbcScatterHighSubtract[number&3];
2359  }
2360#endif
2361  scatter[pivotRow].offset=lastEntryByColumnUPlus_;
2362  scatter[pivotRow].number=number;
2363#endif
2364  // if scatter could move from inside pivot etc rather than indexRow
2365  CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
2366  CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
2367  CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
2368  CoinBigIndex startC = startColumnU[pivotColumn]-number;
2369  CoinFactorizationDouble pivotMultiplier = pivotRegionAddress_[numberGoodU_];
2370  for (int i=startC;i<startC+number;i++)
2371    elementUColumnPlusAddress_[lastEntryByColumnUPlus_++]=element[i]*pivotMultiplier;
2372  CoinAbcMemcpy(reinterpret_cast<CoinSimplexInt *>(elementUColumnPlusAddress_+lastEntryByColumnUPlus_),indexRow+startC,number);
2373  lastEntryByColumnUPlus_ += (number+1)>>1;
2374#ifdef ABC_USE_FUNCTION_POINTERS
2375  if (!numberInColumnAddress_[pivotColumn]) {
2376    int iNext=nextColumnAddress_[pivotColumn];
2377    int iLast=lastColumnAddress_[pivotColumn];
2378    lastColumnAddress_[iNext]=iLast;
2379    assert (iLast>=0);
2380    nextColumnAddress_[iLast]=iNext;
2381  }
2382#endif
2383#endif
2384  numberGoodU_++;
2385  numberRowsLeft_--;
2386#if COIN_DEBUG==2
2387  checkConsistency (  );
2388#endif
2389}
2390#if ABC_DENSE_CODE
2391// After pivoting - returns true if need to go dense
2392int CoinAbcTypeFactorization::wantToGoDense()
2393{
2394  int status=0;
2395  if (denseThreshold_) {
2396    // see whether to go dense
2397    // only if exact multiple
2398    if ((numberRowsLeft_&7)==0) {
2399      CoinSimplexDouble full = numberRowsLeft_;
2400      full *= full;
2401      assert (full>=0.0);
2402      CoinSimplexDouble leftElements = totalElements_;
2403      //if (numberRowsLeft_==100)
2404      //printf("at 100 %d elements\n",totalElements_);
2405      CoinSimplexDouble ratio;
2406      if (numberRowsLeft_>2000)
2407        ratio=4.5;
2408      else if (numberRowsLeft_>800)
2409        ratio=3.0;//3.5;
2410      else if (numberRowsLeft_>256)
2411        ratio=2.0;//2.5;
2412      //else if (numberRowsLeft_==8)
2413      //ratio=200.0;
2414      else
2415        ratio=1.5;
2416      //ratio *=0.75;
2417      if ((ratio*leftElements>full&&numberRowsLeft_>abs(denseThreshold_))) {
2418#if 0 //ndef NDEBUG
2419        //return to do dense
2420        CoinSimplexInt check=0;
2421        for (CoinSimplexInt iColumn=0;iColumn<numberRowsSmall_;iColumn++) {
2422          if (numberInColumnAddress_[iColumn])
2423            check++;
2424        }
2425        if (check!=numberRowsLeft_) {
2426          printf("** mismatch %d columns left, %d rows\n",check,numberRowsLeft_);
2427          //denseThreshold_=0;
2428        } else {
2429          status=2;
2430          if ((messageLevel_&4)!=0)
2431            std::cout<<"      Went dense at "<<numberRowsLeft_<<" rows "<<
2432              totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
2433        }
2434#else
2435        status=2;
2436#if 0 //ABC_NORMAL_DEBUG>1
2437        std::cout<<"      Went dense at "<<numberRowsLeft_<<" rows "<<
2438          totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
2439#endif
2440#endif
2441      }
2442    }
2443  }
2444  return status;
2445}
2446#endif
2447//  pivotRowSingleton.  Does one pivot on Row Singleton in factorization
2448bool
2449CoinAbcTypeFactorization::pivotRowSingleton ( CoinSimplexInt pivotRow,
2450                                  CoinSimplexInt pivotColumn )
2451{
2452  //store pivot columns (so can easily compress)
2453  CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
2454  CoinBigIndex  COIN_RESTRICT startColumn = startColumnU[pivotColumn];
2455  CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
2456  CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
2457  CoinSimplexInt numberDoColumn = numberInColumn[pivotColumn] - 1;
2458  CoinBigIndex endColumn = startColumn + numberDoColumn + 1;
2459  CoinBigIndex pivotRowPosition = startColumn;
2460  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
2461  CoinSimplexInt iRow = indexRowU[pivotRowPosition];
2462  CoinBigIndex *  COIN_RESTRICT startRowU = startRowUAddress_;
2463  CoinSimplexInt *  COIN_RESTRICT nextRow = nextRowAddress_;
2464  CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
2465
2466  while ( iRow != pivotRow ) {
2467    pivotRowPosition++;
2468    iRow = indexRowU[pivotRowPosition];
2469  }                             /* endwhile */
2470  assert ( pivotRowPosition < endColumn );
2471  //store column in L, compress in U and take column out
2472  CoinBigIndex l = lengthL_;
2473
2474  if ( l + numberDoColumn > lengthAreaL_ ) {
2475    //need more memory
2476    if ((messageLevel_&4)!=0) 
2477      std::cout << "more memory needed in middle of invert" << std::endl;
2478    return false;
2479  }
2480  CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
2481  CoinFactorizationDouble *  COIN_RESTRICT elementL = elementLAddress_;
2482  CoinSimplexInt *  COIN_RESTRICT indexRowL = indexRowLAddress_;
2483  startColumnL[numberGoodL_] = l;       //for luck and first time
2484  numberGoodL_++;
2485  startColumnL[numberGoodL_] = l + numberDoColumn;
2486  lengthL_ += numberDoColumn;
2487  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
2488  CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
2489  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
2490
2491  pivotRegionAddress_[numberGoodU_] = pivotMultiplier;
2492  CoinBigIndex i;
2493
2494  CoinSimplexInt *  COIN_RESTRICT indexColumnU = indexColumnUAddress_;
2495#ifdef SMALL_PERMUTE
2496  int realPivotRow=fromSmallToBigRow_[pivotRow];
2497#endif
2498  for ( i = startColumn; i < pivotRowPosition; i++ ) {
2499    CoinSimplexInt iRow = indexRowU[i];
2500
2501#ifdef SMALL_PERMUTE
2502    indexRowL[l] = fromSmallToBigRow_[iRow];
2503#else
2504    indexRowL[l] = iRow;
2505#endif
2506    elementL[l] = elementU[i] * pivotMultiplier;
2507    l++;
2508    //take out of row list
2509    CoinBigIndex start = startRowU[iRow];
2510    CoinSimplexInt iNumberInRow = numberInRow[iRow];
2511    CoinBigIndex end = start + iNumberInRow;
2512    CoinBigIndex where = start;
2513
2514    while ( indexColumnU[where] != pivotColumn ) {
2515      where++;
2516    }                           /* endwhile */
2517    assert ( where < end );
2518#if BOTH_WAYS
2519#endif
2520    indexColumnU[where] = indexColumnU[end - 1];
2521    iNumberInRow--;
2522    numberInRow[iRow] = iNumberInRow;
2523    modifyLink ( iRow, iNumberInRow );
2524  }
2525  for ( i = pivotRowPosition + 1; i < endColumn; i++ ) {
2526    CoinSimplexInt iRow = indexRowU[i];
2527
2528#ifdef SMALL_PERMUTE
2529    indexRowL[l] = fromSmallToBigRow_[iRow];
2530#else
2531    indexRowL[l] = iRow;
2532#endif
2533    elementL[l] = elementU[i] * pivotMultiplier;
2534    l++;
2535    //take out of row list
2536    CoinBigIndex start = startRowU[iRow];
2537    CoinSimplexInt iNumberInRow = numberInRow[iRow];
2538    CoinBigIndex end = start + iNumberInRow;
2539    CoinBigIndex where = start;
2540
2541    while ( indexColumnU[where] != pivotColumn ) {
2542      where++;
2543    }                           /* endwhile */
2544    assert ( where < end );
2545#if BOTH_WAYS
2546#endif
2547    indexColumnU[where] = indexColumnU[end - 1];
2548    iNumberInRow--;
2549    numberInRow[iRow] = iNumberInRow;
2550    modifyLink ( iRow, iNumberInRow );
2551  }
2552  numberInColumn[pivotColumn] = 0;
2553  //modify linked list for pivots
2554  numberInRow[pivotRow] = 0;
2555  deleteLink ( pivotRow );
2556  deleteLink ( pivotColumn + numberRows_ );
2557  //take out this bit of indexColumnU
2558  CoinSimplexInt next = nextRow[pivotRow];
2559  CoinSimplexInt last = lastRow[pivotRow];
2560
2561  nextRow[last] = next;
2562  lastRow[next] = last;
2563  lastRow[pivotRow] =-2;
2564#ifdef SMALL_PERMUTE
2565  permuteAddress_[realPivotRow]=numberGoodU_;
2566#else
2567  permuteAddress_[pivotRow]=numberGoodU_;
2568#endif
2569  return true;
2570}
2571
2572//  pivotColumnSingleton.  Does one pivot on Column Singleton in factorization
2573void
2574CoinAbcTypeFactorization::pivotColumnSingleton ( CoinSimplexInt pivotRow,
2575                                     CoinSimplexInt pivotColumn )
2576{
2577  CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
2578  CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
2579  CoinSimplexInt *  COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
2580  //store pivot columns (so can easily compress)
2581  CoinSimplexInt numberDoRow = numberInRow[pivotRow] - 1;
2582  CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
2583  CoinBigIndex startColumn = startColumnU[pivotColumn];
2584  CoinBigIndex *  COIN_RESTRICT startRowU = startRowUAddress_;
2585  CoinBigIndex startRow = startRowU[pivotRow];
2586  CoinBigIndex endRow = startRow + numberDoRow + 1;
2587  CoinSimplexInt *  COIN_RESTRICT nextRow = nextRowAddress_;
2588  CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
2589  //take out this bit of indexColumnU
2590  CoinSimplexInt next = nextRow[pivotRow];
2591  CoinSimplexInt last = lastRow[pivotRow];
2592
2593  nextRow[last] = next;
2594  lastRow[next] = last;
2595  lastRow[pivotRow] =-2; //mark
2596#ifdef SMALL_PERMUTE
2597  int realPivotRow=fromSmallToBigRow_[pivotRow];
2598  permuteAddress_[realPivotRow]=numberGoodU_;
2599#else
2600  permuteAddress_[pivotRow]=numberGoodU_;
2601#endif
2602  //clean up counts
2603  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
2604  CoinFactorizationDouble pivotElement = elementU[startColumn];
2605
2606  pivotRegionAddress_[numberGoodU_] = 1.0 / pivotElement;
2607  numberInColumn[pivotColumn] = 0;
2608  CoinSimplexInt *  COIN_RESTRICT indexColumnU = indexColumnUAddress_;
2609  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
2610  //CoinSimplexInt *  COIN_RESTRICT saveColumn = saveColumnAddress_;
2611  CoinBigIndex i;
2612  // would I be better off doing other way first??
2613  for ( i = startRow; i < endRow; i++ ) {
2614    CoinSimplexInt iColumn = indexColumnU[i];
2615    if ( numberInColumn[iColumn] ) {
2616      CoinSimplexInt number = numberInColumn[iColumn] - 1;
2617      //modify linked list
2618      modifyLink ( iColumn + numberRows_, number );
2619      //move pivot row element
2620#ifdef SMALL_PERMUTE
2621      CoinBigIndex start = startColumnU[iColumn];
2622#endif
2623      if ( number ) {
2624#ifndef SMALL_PERMUTE
2625        CoinBigIndex start = startColumnU[iColumn];
2626#endif
2627        CoinBigIndex pivot = start;
2628        CoinSimplexInt iRow = indexRowU[pivot];
2629        while ( iRow != pivotRow ) {
2630          pivot++;
2631          iRow = indexRowU[pivot];
2632        }
2633        assert (pivot < startColumnU[iColumn] +
2634                numberInColumn[iColumn]);
2635#if BOTH_WAYS
2636#endif
2637        if ( pivot != start ) {
2638          //move largest one up
2639          CoinFactorizationDouble value = elementU[start];
2640
2641          iRow = indexRowU[start];
2642          elementU[start] = elementU[pivot];
2643          assert (indexRowU[pivot]==pivotRow);
2644          indexRowU[start] = indexRowU[pivot];
2645          elementU[pivot] = elementU[start + 1];
2646          indexRowU[pivot] = indexRowU[start + 1];
2647          elementU[start + 1] = value;
2648          indexRowU[start + 1] = iRow;
2649#if BOTH_WAYS
2650#endif
2651        } else {
2652          //find new largest element
2653          CoinSimplexInt iRowSave = indexRowU[start + 1];
2654          CoinFactorizationDouble valueSave = elementU[start + 1];
2655          CoinFactorizationDouble valueLargest = fabs ( valueSave );
2656          CoinBigIndex end = start + numberInColumn[iColumn];
2657          CoinBigIndex largest = start + 1;
2658
2659          CoinBigIndex k;
2660          for ( k = start + 2; k < end; k++ ) {
2661            CoinFactorizationDouble value = elementU[k];
2662            CoinFactorizationDouble valueAbs = fabs ( value );
2663
2664            if ( valueAbs > valueLargest ) {
2665              valueLargest = valueAbs;
2666              largest = k;
2667            }
2668          }
2669          indexRowU[start + 1] = indexRowU[largest];
2670          elementU[start + 1] = elementU[largest];
2671          indexRowU[largest] = iRowSave;
2672          elementU[largest] = valueSave;
2673#if BOTH_WAYS
2674#endif
2675        }
2676      }
2677#ifdef SMALL_PERMUTE
2678      indexRowU[start]=realPivotRow;
2679#endif
2680      //clean up counts
2681      numberInColumn[iColumn]--;
2682      numberInColumnPlus[iColumn]++;
2683      startColumnU[iColumn]++;
2684      //totalElements_--;
2685    }
2686  }
2687  //modify linked list for pivots
2688  deleteLink ( pivotRow );
2689  deleteLink ( pivotColumn + numberRows_ );
2690  numberInRow[pivotRow] = 0;
2691  //put in dummy pivot in L
2692  CoinBigIndex l = lengthL_;
2693
2694  CoinBigIndex * startColumnL = startColumnLAddress_;
2695  startColumnL[numberGoodL_] = l;       //for luck and first time
2696  numberGoodL_++;
2697  startColumnL[numberGoodL_] = l;
2698}
2699
2700
2701//  getColumnSpace.  Gets space for one Column with given length
2702//may have to do compression  (returns true)
2703//also moves existing vector
2704bool
2705CoinAbcTypeFactorization::getColumnSpace ( CoinSimplexInt iColumn,
2706                               CoinSimplexInt extraNeeded )
2707{
2708  CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
2709  CoinSimplexInt *  COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
2710  CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumnAddress_;
2711  CoinSimplexInt *  COIN_RESTRICT lastColumn = lastColumnAddress_;
2712  CoinSimplexInt number = numberInColumnPlus[iColumn] +
2713    numberInColumn[iColumn];
2714  CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
2715  CoinBigIndex space = lengthAreaU_ - lastEntryByColumnU_;
2716  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
2717  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
2718
2719  if ( space < extraNeeded + number + 4 ) {
2720    //compression
2721    CoinSimplexInt iColumn = nextColumn[maximumRowsExtra_];
2722    CoinBigIndex put = 0;
2723
2724    while ( iColumn != maximumRowsExtra_ ) {
2725      //move
2726      CoinBigIndex get;
2727      CoinBigIndex getEnd;
2728
2729      if ( startColumnU[iColumn] >= 0 ) {
2730        get = startColumnU[iColumn]
2731          - numberInColumnPlus[iColumn];
2732        getEnd = startColumnU[iColumn] + numberInColumn[iColumn];
2733        startColumnU[iColumn] = put + numberInColumnPlus[iColumn];
2734      } else {
2735        get = -startColumnU[iColumn];
2736        getEnd = get + numberInColumn[iColumn];
2737        startColumnU[iColumn] = -put;
2738      }
2739      assert (put<=get);
2740      for (CoinBigIndex i = get; i < getEnd; i++ ) {
2741        indexRowU[put] = indexRowU[i];
2742        elementU[put] = elementU[i];
2743#if BOTH_WAYS
2744#endif
2745        put++;
2746      }
2747      iColumn = nextColumn[iColumn];
2748    }                           /* endwhile */
2749    numberCompressions_++;
2750    lastEntryByColumnU_ = put;
2751    space = lengthAreaU_ - put;
2752    if ( extraNeeded == COIN_INT_MAX >> 1 ) {
2753      return true;
2754    }
2755    if ( space < extraNeeded + number + 2 ) {
2756      //need more space
2757      //if we can allocate bigger then do so and copy
2758      //if not then return so code can start again
2759      status_ = -99;
2760      return false;
2761    }
2762  }
2763  CoinBigIndex put = lastEntryByColumnU_;
2764  CoinSimplexInt next = nextColumn[iColumn];
2765  CoinSimplexInt last = lastColumn[iColumn];
2766
2767  if ( extraNeeded || next != maximumRowsExtra_ ) {
2768    //out
2769    nextColumn[last] = next;
2770    lastColumn[next] = last;
2771    //in at end
2772    last = lastColumn[maximumRowsExtra_];
2773    nextColumn[last] = iColumn;
2774    lastColumn[maximumRowsExtra_] = iColumn;
2775    lastColumn[iColumn] = last;
2776    nextColumn[iColumn] = maximumRowsExtra_;
2777    //move
2778    CoinBigIndex get = startColumnU[iColumn]
2779      - numberInColumnPlus[iColumn];
2780
2781    startColumnU[iColumn] = put + numberInColumnPlus[iColumn];
2782    CoinSimplexInt *indexRow = indexRowU;
2783    CoinFactorizationDouble *element = elementU;
2784    CoinSimplexInt i = 0;
2785   
2786    if ( ( number & 1 ) != 0 ) {
2787      element[put] = element[get];
2788      indexRow[put] = indexRow[get];
2789#if BOTH_WAYS
2790#endif
2791      i = 1;
2792    }
2793    for ( ; i < number; i += 2 ) {
2794      CoinFactorizationDouble value0 = element[get + i];
2795      CoinFactorizationDouble value1 = element[get + i + 1];
2796      CoinSimplexInt index0 = indexRow[get + i];
2797      CoinSimplexInt index1 = indexRow[get + i + 1];
2798     
2799      element[put + i] = value0;
2800      element[put + i + 1] = value1;
2801      indexRow[put + i] = index0;
2802      indexRow[put + i + 1] = index1;
2803#if BOTH_WAYS
2804#endif
2805    }
2806    put += number;
2807    get += number;
2808    //add 2 for luck
2809    lastEntryByColumnU_ = put + extraNeeded + 2;
2810    if (lastEntryByColumnU_>lengthAreaU_) {
2811      // get more memory
2812#ifdef CLP_DEVELOP
2813      printf("put %d, needed %d, start %d, length %d\n",
2814             put,extraNeeded,lastEntryByColumnU_,
2815             lengthAreaU_);
2816#endif
2817      return false;
2818    }
2819  } else {
2820    //take off space
2821    lastEntryByColumnU_ = startColumnU[last] + numberInColumn[last];
2822  }
2823  startColumnU[maximumRowsExtra_]=lastEntryByColumnU_;
2824  return true;
2825}
2826
2827//  getRowSpace.  Gets space for one Row with given length
2828//may have to do compression  (returns true)
2829//also moves existing vector
2830bool
2831CoinAbcTypeFactorization::getRowSpace ( CoinSimplexInt iRow,
2832                            CoinSimplexInt extraNeeded )
2833{
2834  CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
2835  CoinSimplexInt number = numberInRow[iRow];
2836  CoinBigIndex *  COIN_RESTRICT startRowU = startRowUAddress_;
2837  CoinBigIndex space = lengthAreaU_ - lastEntryByRowU_;
2838  CoinSimplexInt *  COIN_RESTRICT nextRow = nextRowAddress_;
2839  CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
2840  CoinSimplexInt *  COIN_RESTRICT indexColumnU = indexColumnUAddress_;
2841
2842  if ( space < extraNeeded + number + 2 ) {
2843    //compression
2844    CoinSimplexInt iRow = nextRow[numberRows_];
2845    CoinBigIndex put = 0;
2846
2847    while ( iRow != numberRows_ ) {
2848      //move
2849      CoinBigIndex get = startRowU[iRow];
2850      CoinBigIndex getEnd = startRowU[iRow] + numberInRow[iRow];
2851
2852      startRowU[iRow] = put;
2853      if (put>get) {
2854        //need more space
2855        //if we can allocate bigger then do so and copy
2856        //if not then return so code can start again
2857        status_ = -99;
2858        return false;;
2859      }
2860      CoinBigIndex i;
2861      for ( i = get; i < getEnd; i++ ) {
2862        indexColumnU[put] = indexColumnU[i];
2863#if BOTH_WAYS
2864#endif
2865        put++;
2866      }
2867      iRow = nextRow[iRow];
2868    }                           /* endwhile */
2869    numberCompressions_++;
2870    lastEntryByRowU_ = put;
2871    space = lengthAreaU_ - put;
2872    if ( space < extraNeeded + number + 2 ) {
2873      //need more space
2874      //if we can allocate bigger then do so and copy
2875      //if not then return so code can start again
2876      status_ = -99;
2877      return false;;
2878    }
2879  }
2880  CoinBigIndex put = lastEntryByRowU_;
2881  CoinSimplexInt next = nextRow[iRow];
2882  CoinSimplexInt last = lastRow[iRow];
2883
2884  //out
2885  nextRow[last] = next;
2886  lastRow[next] = last;
2887  //in at end
2888  last = lastRow[numberRows_];
2889  nextRow[last] = iRow;
2890  lastRow[numberRows_] = iRow;
2891  lastRow[iRow] = last;
2892  nextRow[iRow] = numberRows_;
2893  //move
2894  CoinBigIndex get = startRowU[iRow];
2895
2896  startRowU[iRow] = put;
2897  while ( number ) {
2898    number--;
2899    indexColumnU[put] = indexColumnU[get];
2900#if BOTH_WAYS
2901#endif
2902    put++;
2903    get++;
2904  }                             /* endwhile */
2905  //add 4 for luck
2906  lastEntryByRowU_ = put + extraNeeded + 4;
2907  return true;
2908}
2909
2910//  cleanup.  End of factorization
2911void
2912CoinAbcTypeFactorization::cleanup (  )
2913{
2914#ifdef ABC_USE_FUNCTION_POINTERS
2915#define PRINT_LEVEL 0
2916  if (PRINT_LEVEL){
2917    printf("----------Start---------\n");
2918#if PRINT_LEVEL<2
2919    double smallestP=COIN_DBL_MAX;
2920    int iPs=-1;
2921    double largestP=0.0;
2922    int iP=-1;
2923    double largestL=0.0;
2924    int iL=-1;
2925    double largestU=0.0;
2926    int iU=-1;
2927#endif
2928    for (int i=0;i<numberRows_;i++) {
2929#if PRINT_LEVEL>1
2930      printf("%d permute %d pivotColumn %d pivotRegion %g\n",i,permuteAddress_[i],
2931             pivotColumnAddress_[i],pivotRegionAddress_[i]);
2932#else
2933      if (fabs(pivotRegionAddress_[i])>largestP) {
2934        iP=i;
2935        largestP=fabs(pivotRegionAddress_[i]);
2936      }
2937      if (fabs(pivotRegionAddress_[i])<smallestP) {
2938        iPs=i;
2939        smallestP=fabs(pivotRegionAddress_[i]);
2940      }
2941#endif
2942    }
2943    scatterStruct * COIN_RESTRICT scatterPointer = scatterUColumn();
2944    CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
2945    for (int i=0;i<numberRows_;i++) {
2946      scatterStruct & COIN_RESTRICT scatter = scatterPointer[i];
2947      CoinSimplexInt number = scatter.number;
2948#if PRINT_LEVEL>1
2949      printf("pivot %d offset %d number %d\n",i,scatter.offset,number);
2950#endif
2951      const CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter.offset;
2952      const CoinBigIndex * COIN_RESTRICT indices = reinterpret_cast<const CoinBigIndex *>(area+number);
2953#if PRINT_LEVEL>1
2954      for (int j1 = 0; j1 < number; j1+=5 ) {
2955        for (int j = j1; j < CoinMin(number,j1+5); j++ ) 
2956          printf("(%d,%g) ",indices[j],area[j]);
2957        printf("\n");
2958      }
2959#else
2960      for (int j = 0; j < number; j++ ) 
2961        if (fabs(area[j])>largestU) {
2962          iU=i;
2963          largestU=fabs(area[j]);
2964        }
2965#endif
2966    }
2967    printf("----------LLL---------\n");
2968    CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
2969    CoinSimplexInt *  COIN_RESTRICT indexRowL = indexRowLAddress_;
2970    CoinFactorizationDouble *  COIN_RESTRICT elementL = elementLAddress_;
2971    for (int i=0;i<numberRows_;i++) {
2972      int number=startColumnL[i+1]-startColumnL[i];
2973      if (number) {
2974#if PRINT_LEVEL>1
2975        int * indices = indexRowL+startColumnL[i];
2976        printf("%d has %d elements\n",i,number);
2977        for (int j1 = 0; j1 < number; j1+=5 ) {
2978          for (int j = j1; j < CoinMin(number,j1+5); j++ ) 
2979            printf("(%d,%g) ",indices[j],elementL[j+startColumnL[i]]);
2980          printf("\n");
2981        }
2982#else
2983        for (int j = 0; j < number; j++ ) {
2984          if (fabs(elementL[j+startColumnL[i]])>largestL) {
2985            iL=i;
2986            largestL=fabs(elementL[j+startColumnL[i]]);
2987          }
2988        }
2989#endif
2990      }
2991    }
2992#if PRINT_LEVEL<2
2993    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);
2994#endif
2995    printf("----------End---------\n");
2996  }
2997#endif
2998#ifdef ABC_ORDERED_FACTORIZATION
2999  {
3000    //int sizeL=0;
3001    //int sizeU=0;
3002    CoinSimplexInt *  COIN_RESTRICT permuteA = permuteAddress_+maximumRowsExtra_+1;
3003    CoinSimplexInt *  COIN_RESTRICT permuteB = permuteA + numberRows_;
3004    // one may be correct
3005    for (int i=0;i<numberRows_;i++) {
3006      int iPermute=permuteAddress_[i];
3007      permuteA[i]=iPermute;
3008      permuteB[iPermute]=i;
3009      permuteAddress_[i]=i;
3010    }
3011    scatterStruct * COIN_RESTRICT scatterPointer = scatterUColumn();
3012    CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
3013    for (int i=0;i<numberRows_;i++) {
3014      scatterStruct & COIN_RESTRICT scatter = scatterPointer[i];
3015      CoinSimplexInt number = scatter.number;
3016      //sizeU+=number;
3017      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter.offset;
3018      CoinBigIndex * COIN_RESTRICT indices = reinterpret_cast<CoinBigIndex *>(area+number);
3019      for (int j = 0; j < number; j++ ) {
3020        int iRow=indices[j];
3021        indices[j]=permuteA[iRow];
3022      }
3023    }
3024    CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
3025    CoinSimplexInt *  COIN_RESTRICT indexRowL = indexRowLAddress_;
3026    CoinFactorizationDouble *  COIN_RESTRICT elementL = elementLAddress_;
3027    for (int i=0;i<numberRows_;i++) {
3028      int number=startColumnL[i+1]-startColumnL[i];
3029      if (number) {
3030        //sizeL+=number;
3031        int * indices = indexRowL+startColumnL[i];
3032        for (int j = 0; j < number; j++ ) {
3033          int iRow=indices[j];
3034          indices[j]=permuteA[iRow];
3035        }
3036      }
3037    }
3038    //printf("Ordered says %d L %d U\n",sizeL,sizeU);
3039  }
3040#elif 0
3041#ifdef ABC_USE_FUNCTION_POINTERS
3042  {
3043    int sizeL=0;
3044    int sizeU=0;
3045    scatterStruct * COIN_RESTRICT scatterPointer = scatterUColumn();
3046    for (int i=0;i<numberRows_;i++) {
3047      scatterStruct & COIN_RESTRICT scatter = scatterPointer[i];
3048      CoinSimplexInt number = scatter.number;
3049      sizeU+=number;
3050    }
3051    CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
3052    for (int i=0;i<numberRows_;i++) {
3053      int number=startColumnL[i+1]-startColumnL[i];
3054      if (number) 
3055        sizeL+=number;
3056    }
3057    printf("Unordered says %d L %d U\n",sizeL,sizeU);
3058  }
3059#endif
3060#endif
3061  //free some memory here
3062  markRow_.conditionalDelete() ;
3063  nextRowAddress_ = nextRow_.array();
3064  //safety feature
3065  CoinSimplexInt *  COIN_RESTRICT permute = permuteAddress_;
3066  permute[numberRows_] = 0;
3067#ifndef NDEBUG
3068  {
3069    if (numberGoodU_<numberRows_)
3070      abort();
3071    char * mark = new char[numberRows_];
3072    memset(mark,0,numberRows_);
3073    for (int i = 0; i < numberRows_; i++ ) {
3074      assert (permute[i]>=0&&permute[i]<numberRows_);
3075      CoinSimplexInt k = permute[i];
3076      if(k<0||k>=numberRows_)
3077        printf("Bad a %d %d\n",i,k);
3078      assert(k>=0&&k<numberRows_);
3079      if(mark[k]==1)
3080        printf("Bad a %d %d\n",i,k);
3081      mark[k]=1;
3082    }
3083    for (int i = 0; i < numberRows_; i++ ) {
3084      assert(mark[i]==1);
3085      if(mark[i]!=1)
3086        printf("Bad b %d\n",i);
3087    }
3088    delete [] mark;
3089  }
3090#endif
3091  CoinSimplexInt *  COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
3092  const CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_;
3093  CoinSimplexInt *  COIN_RESTRICT pivotLOrder = firstCountAddress_;//tttt;//this->pivotLOrder();
3094  CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
3095  bool cleanCopy=false;
3096#if ABC_SMALL<2
3097  assert (numberGoodU_==numberRows_);
3098  if (numberGoodU_==numberRows_) {
3099#ifndef ABC_USE_FUNCTION_POINTERS
3100    const CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnPlusAddress_;
3101    CoinFactorizationDouble * elementU = elementUAddress_;
3102    CoinFactorizationDouble * elementU2 = elementRowUAddress_;
3103    CoinSimplexInt * COIN_RESTRICT indexRowU = indexRowUAddress_; 
3104    CoinSimplexInt * COIN_RESTRICT indexRowU2 = indexColumnUAddress_; 
3105    // so that just slacks will have start of 0
3106    CoinBigIndex size=1;
3107    for (int i=0;i<numberSlacks_;i++) {
3108      CoinSimplexInt iNext=pivotColumn[i];
3109      startColumnU[iNext]=0;
3110      assert (!numberInColumn[iNext]);
3111    }
3112    for (int i=numberSlacks_;i<numberRows_;i++) {
3113      CoinSimplexInt iNext=pivotColumn[i];
3114      CoinSimplexInt number = numberInColumn[iNext];
3115      CoinBigIndex start = startColumnU[iNext]-number;
3116      startColumnU[iNext]=size;
3117      CoinAbcMemcpy(elementU2+size,elementU+start,number);
3118      CoinAbcMemcpy(indexRowU2+size,indexRowU+start,number);
3119      size+=number;
3120    }
3121    // swap arrays
3122    elementU_.swap(elementRowU_);
3123    elementUAddress_ = elementU_.array();
3124    elementRowUAddress_ = elementRowU_.array();
3125    indexColumnU_.swap(indexRowU_);
3126    indexRowUAddress_ = indexRowU_.array();
3127    indexColumnUAddress_ = indexColumnU_.array();
3128    // Redo total elements
3129    totalElements_=size;
3130#else
3131    // swap arrays
3132    elementU_.swap(elementRowU_);
3133    elementUAddress_ = elementU_.array();
3134    elementRowUAddress_ = elementRowU_.array();
3135    totalElements_=lastEntryByColumnUPlus_;
3136#endif
3137    cleanCopy=true;
3138  }
3139#endif
3140#ifndef ABC_USE_FUNCTION_POINTERS
3141  if (!cleanCopy) {
3142    getColumnSpace ( 0, COIN_INT_MAX >> 1 );    //compress
3143    // Redo total elements
3144    totalElements_=0;
3145    for (int i = 0; i < numberRows_; i++ ) {
3146      CoinSimplexInt number = numberInColumnPlus[i];
3147      totalElements_ += number;
3148      startColumnU[i] -= number;
3149    }
3150  }
3151#endif
3152  // swap arrays
3153  numberInColumn_.swap(numberInColumnPlus_);
3154  numberInColumnAddress_ = numberInColumn_.array();
3155  numberInColumnPlusAddress_ = numberInColumnPlus_.array();
3156  CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
3157  CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
3158#if ABC_SMALL<2
3159  numberInColumnPlus = numberInColumnPlusAddress_;
3160#endif
3161  //make column starts OK
3162  //for best cache behavior get in order (last pivot at bottom of space)
3163  //that will need thinking about
3164
3165  lastEntryByRowU_ = totalElements_;
3166
3167  CoinSimplexInt *  COIN_RESTRICT back = firstCountAddress_+numberRows_;
3168  CoinBigIndex * COIN_RESTRICT startRow = startRowUAddress_;
3169#ifndef ABC_USE_FUNCTION_POINTERS
3170  maximumU_=startColumnU[numberRows_-1]+numberInColumn[numberRows_-1];
3171#else
3172  maximumU_=lastEntryByColumnUPlus_;
3173#endif
3174#ifndef ABC_ORDERED_FACTORIZATION
3175  CoinFactorizationDouble *  COIN_RESTRICT workArea = workAreaAddress_;
3176  CoinFactorizationDouble *  COIN_RESTRICT pivotRegion = pivotRegionAddress_;
3177#else
3178  CoinSimplexInt *  COIN_RESTRICT permuteA = permuteAddress_+maximumRowsExtra_+1;
3179  // use work area for scatter
3180  scatterStruct * COIN_RESTRICT scatterPointer = scatterUColumn();
3181  scatterStruct * COIN_RESTRICT scatterPointerTemp = reinterpret_cast<scatterStruct *>(workAreaAddress_);
3182#endif
3183  for (CoinSimplexInt i=0;i<numberRows_;i++) {
3184#ifndef ABC_ORDERED_FACTORIZATION
3185    CoinSimplexInt pivotOrder=permute[i];
3186    CoinSimplexInt columnOrder=pivotColumn[pivotOrder];
3187    workArea[i]=pivotRegion[pivotOrder];
3188    back[columnOrder] = i;
3189#else
3190    CoinSimplexInt pivotOrder=i;
3191    CoinSimplexInt columnOrder=pivotColumn[pivotOrder];
3192    int pivotOrder2=permuteA[i];
3193    scatterPointerTemp[pivotOrder2]=scatterPointer[i];
3194    back[columnOrder] = i;
3195#endif
3196#ifndef ABC_USE_FUNCTION_POINTERS
3197    startRow[i]=startColumnU[columnOrder];
3198    numberInRow[i]=numberInColumn[columnOrder];
3199#endif
3200    pivotLOrder[pivotOrder]=i;
3201  }
3202#ifndef ABC_ORDERED_FACTORIZATION
3203  CoinAbcMemcpy(pivotRegion,workArea,numberRows_); // could swap
3204#else
3205  CoinAbcMemcpy(scatterPointer,scatterPointerTemp,numberRows_);
3206#endif
3207#ifndef EARLY_FACTORIZE
3208  workArea_.conditionalDelete() ;
3209  workAreaAddress_=NULL;
3210#endif
3211#ifndef ABC_USE_FUNCTION_POINTERS
3212  //memcpy(this->pivotLOrder(),pivotLOrder,numberRows_*sizeof(int));
3213  // swap later
3214  CoinAbcMemcpy(startColumnU,startRow,numberRows_);
3215  CoinAbcMemcpy(numberInColumn,numberInRow,numberRows_);
3216#endif
3217#if ABC_DENSE_CODE==2
3218  if (numberDense_) {
3219    assert (numberDense_<30000);
3220    CoinFactorizationDouble *  COIN_RESTRICT denseArea = denseAreaAddress_;
3221    CoinFactorizationDouble * COIN_RESTRICT denseRegion = 
3222      denseArea+leadingDimension_*numberDense_;
3223    CoinSimplexInt *  COIN_RESTRICT densePermute=
3224      reinterpret_cast<CoinSimplexInt *>(denseRegion+FACTOR_CPU*numberDense_);
3225    short * COIN_RESTRICT forFtran = reinterpret_cast<short *>(densePermute+numberDense_);
3226    short * COIN_RESTRICT forBtran = forFtran+numberDense_;
3227    short * COIN_RESTRICT forFtran2 = forBtran+numberDense_;
3228    // maybe could do inside dgetrf
3229    //CoinAbcMemcpy(forBtran,pivotLOrder+numberRows_-numberDense_,numberDense_);
3230    //CoinAbcMemcpy(forFtran,pivotLOrder+numberRows_-numberDense_,numberDense_);
3231    for (int i=0;i<numberDense_;i++) {
3232      forBtran[i]=static_cast<short>(i);
3233      forFtran[i]=static_cast<short>(i);
3234    }
3235    for (int i=0;i<numberDense_;i++) {
3236      int ip=densePermute[i];
3237      short temp=forBtran[i];
3238      forBtran[i]=forBtran[ip];
3239      forBtran[ip]=temp;
3240    }
3241    for (int i=numberDense_-1;i>=0;i--) {
3242      int ip=densePermute[i];
3243      short temp=forFtran[i];
3244      forFtran[i]=forFtran[ip];
3245      forFtran[ip]=temp;
3246    }
3247#if ABC_SMALL<3
3248    const CoinSimplexInt * COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
3249    int lastSparse=numberRows_-numberDense_;
3250    forFtran -= lastSparse; // adjust
3251    CoinFillN(forFtran2,numberRows_,static_cast<short>(-1));
3252    for (int i=0;i<numberRows_;i++) {
3253      int jRow=pivotLBackwardOrder[i];
3254      if (jRow>=lastSparse) {
3255        short kRow = static_cast<short>(forFtran[jRow]);
3256        forFtran2[i]=kRow;
3257      }
3258    }
3259#endif
3260  }
3261#endif
3262  CoinAbcMemset0(numberInRow,numberRows_+1);
3263  if ( (messageLevel_ & 8)) {
3264    std::cout<<"        length of U "<<totalElements_<<", length of L "<<lengthL_;
3265#if ABC_SMALL<4
3266    if (numberDense_)
3267      std::cout<<" plus "<<numberDense_*numberDense_<<" from "<<numberDense_<<" dense rows";
3268#endif
3269    std::cout<<std::endl;
3270  }
3271#if ABC_SMALL<4
3272  // and add L and dense
3273  totalElements_ += numberDense_*numberDense_+lengthL_;
3274#endif
3275#if ABC_SMALL<2
3276  CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumnAddress_;
3277  CoinSimplexInt *  COIN_RESTRICT lastColumn = lastColumnAddress_;
3278  //goSparse2();
3279  // See whether to have extra copy of R
3280#if ABC_SMALL>=0
3281  if (maximumU_>10*numberRows_||numberRows_<200) {
3282    // NO
3283    numberInColumnPlus_.conditionalDelete() ;
3284    numberInColumnPlusAddress_=NULL;
3285    setNoGotRCopy();
3286  } else {
3287#endif
3288    setYesGotRCopy();
3289    for (int i = 0; i < numberRows_; i++ ) {
3290      lastColumn[i] = i - 1;
3291      nextColumn[i] = i + 1;
3292      numberInColumnPlus[i]=0;
3293#if ABC_SMALL>=0
3294    }
3295#endif
3296    nextColumn[numberRows_ - 1] = maximumRowsExtra_;
3297    lastColumn[maximumRowsExtra_] = numberRows_ - 1;
3298    nextColumn[maximumRowsExtra_] = 0;
3299    lastColumn[0] = maximumRowsExtra_;
3300  }
3301#endif
3302  numberU_ = numberGoodU_;
3303  numberL_ = numberGoodL_;
3304
3305  CoinSimplexInt firstReal = numberSlacks_;
3306
3307  CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
3308  // We should know how many in L
3309  //CoinSimplexInt *  COIN_RESTRICT indexRowL = indexRowLAddress_;
3310  for (firstReal = numberSlacks_; firstReal<numberRows_; firstReal++ ) {
3311    if (startColumnL[firstReal+1])
3312      break;
3313  }
3314  totalElements_ += startColumnL[numberRows_];
3315  baseL_ = firstReal;
3316  numberL_ -= firstReal;
3317  factorElements_ = totalElements_;
3318  //use L for R if room
3319  CoinBigIndex space = lengthAreaL_ - lengthL_;
3320  CoinBigIndex spaceUsed = lengthL_ + lengthU_;
3321
3322  CoinSimplexInt needed = ( spaceUsed + numberRows_ - 1 ) / numberRows_;
3323
3324  needed = needed * 2 * maximumPivots_;
3325  if ( needed < 2 * numberRows_ ) {
3326    needed = 2 * numberRows_;
3327  }
3328  if (gotRCopy()) {
3329    //if (numberInColumnPlusAddress_) {
3330    // Need double the space for R
3331    space = space/2;
3332  }
3333  elementRAddress_ = elementLAddress_ + lengthL_;
3334  indexRowRAddress_ = indexRowLAddress_ + lengthL_;
3335  if ( space >= needed ) {
3336    lengthR_ = 0;
3337    lengthAreaR_ = space;
3338  } else {
3339    lengthR_ = 0;
3340    lengthAreaR_ = space;
3341    if ((messageLevel_&4))
3342      std::cout<<"Factorization may need some increasing area space"
3343               <<std::endl;
3344    if ( areaFactor_ ) {
3345      areaFactor_ *= 1.1;
3346    } else {
3347      areaFactor_ = 1.1;
3348    }
3349  }
3350  numberR_ = 0;
3351}
3352// Returns areaFactor but adjusted for dense
3353double 
3354CoinAbcTypeFactorization::adjustedAreaFactor() const
3355{
3356  double factor = areaFactor_;
3357#if ABC_SMALL<4
3358  if (numberDense_&&areaFactor_>1.0) {
3359    double dense = numberDense_;
3360    dense *= dense;
3361    double withoutDense = totalElements_ - dense +1.0;
3362    factor *= 1.0 +dense/withoutDense;
3363  }
3364#endif
3365  return factor;
3366}
3367
3368//  checkConsistency.  Checks that row and column copies look OK
3369void
3370CoinAbcTypeFactorization::checkConsistency (  )
3371{
3372  bool bad = false;
3373
3374  CoinSimplexInt iRow;
3375  CoinBigIndex * startRowU = startRowUAddress_;
3376  CoinSimplexInt * numberInRow = numberInRowAddress_;
3377  CoinSimplexInt * numberInColumn = numberInColumnAddress_;
3378  CoinSimplexInt * indexColumnU = indexColumnUAddress_;
3379  CoinSimplexInt * indexRowU = indexRowUAddress_;
3380  CoinBigIndex * startColumnU = startColumnUAddress_;
3381  for ( iRow = 0; iRow < numberRows_; iRow++ ) {
3382    if ( numberInRow[iRow] ) {
3383      CoinBigIndex startRow = startRowU[iRow];
3384      CoinBigIndex endRow = startRow + numberInRow[iRow];
3385
3386      CoinBigIndex j;
3387      for ( j = startRow; j < endRow; j++ ) {
3388        CoinSimplexInt iColumn = indexColumnU[j];
3389        CoinBigIndex startColumn = startColumnU[iColumn];
3390        CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
3391        bool found = false;
3392
3393        CoinBigIndex k;
3394        for ( k = startColumn; k < endColumn; k++ ) {
3395          if ( indexRowU[k] == iRow ) {
3396            found = true;
3397            break;
3398          }
3399        }
3400        if ( !found ) {
3401          bad = true;
3402          std::cout << "row " << iRow << " column " << iColumn << " Rows" << std::endl;
3403        }
3404      }
3405    }
3406  }
3407  CoinSimplexInt iColumn;
3408  for ( iColumn = 0; iColumn < numberRows_; iColumn++ ) {
3409    if ( numberInColumn[iColumn] ) {
3410      CoinBigIndex startColumn = startColumnU[iColumn];
3411      CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
3412
3413      CoinBigIndex j;
3414      for ( j = startColumn; j < endColumn; j++ ) {
3415        CoinSimplexInt iRow = indexRowU[j];
3416        CoinBigIndex startRow = startRowU[iRow];
3417        CoinBigIndex endRow = startRow + numberInRow[iRow];
3418        bool found = false;
3419
3420        CoinBigIndex k;
3421        for (  k = startRow; k < endRow; k++ ) {
3422          if ( indexColumnU[k] == iColumn ) {
3423            found = true;
3424            break;
3425          }
3426        }
3427        if ( !found ) {
3428          bad = true;
3429          std::cout << "row " << iRow << " column " << iColumn << " Columns" <<
3430            std::endl;
3431        }
3432      }
3433    }
3434  }
3435  if ( bad ) {
3436    abort (  );
3437  }
3438}
3439  //  pivotOneOtherRow.  When just one other row so faster
3440bool 
3441CoinAbcTypeFactorization::pivotOneOtherRow ( CoinSimplexInt pivotRow,
3442                                           CoinSimplexInt pivotColumn )
3443{
3444  CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
3445  CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
3446  CoinSimplexInt *  COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
3447  CoinSimplexInt  COIN_RESTRICT numberInPivotRow = numberInRow[pivotRow] - 1;
3448  CoinBigIndex *  COIN_RESTRICT startRowU = startRowUAddress_;
3449  CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
3450  CoinBigIndex startColumn = startColumnU[pivotColumn];
3451  CoinBigIndex startRow = startRowU[pivotRow];
3452  CoinBigIndex endRow = startRow + numberInPivotRow + 1;
3453
3454  //take out this bit of indexColumnU
3455  CoinSimplexInt *  COIN_RESTRICT nextRow = nextRowAddress_;
3456  CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
3457  CoinSimplexInt next = nextRow[pivotRow];
3458  CoinSimplexInt last = lastRow[pivotRow];
3459
3460  nextRow[last] = next;
3461  lastRow[next] = last;
3462  lastRow[pivotRow] = -2;
3463#ifdef SMALL_PERMUTE
3464  int realPivotRow=fromSmallToBigRow_[pivotRow];
3465  //int realPivotColumn=fromSmallToBigColumn[pivotColumn];
3466#endif
3467#ifdef SMALL_PERMUTE
3468  permuteAddress_[realPivotRow]=numberGoodU_;
3469#else
3470  permuteAddress_[pivotRow]=numberGoodU_;
3471#endif
3472  numberInRow[pivotRow] = 0;
3473  #if ABC_SMALL<2
3474// temp - switch off marker for valid row/column lookup
3475  sparseThreshold_=-2;
3476#endif
3477  //store column in L, compress in U and take column out
3478  CoinBigIndex l = lengthL_;
3479
3480  if ( l + 1 > lengthAreaL_ ) {
3481    //need more memory
3482    if ((messageLevel_&4)!=0) 
3483      std::cout << "more memory needed in middle of invert" << std::endl;
3484    return false;
3485  }
3486  //l+=currentAreaL_->elementByColumn-elementL_;
3487  //CoinBigIndex lSave=l;
3488  CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
3489  CoinFactorizationDouble *  COIN_RESTRICT elementL = elementLAddress_;
3490  CoinSimplexInt *  COIN_RESTRICT indexRowL = indexRowLAddress_;
3491  startColumnL[numberGoodL_] = l;       //for luck and first time
3492  numberGoodL_++;
3493  startColumnL[numberGoodL_] = l + 1;
3494  lengthL_++;
3495  CoinFactorizationDouble pivotElement;
3496  CoinFactorizationDouble otherMultiplier;
3497  CoinSimplexInt otherRow;
3498  CoinSimplexInt *  COIN_RESTRICT saveColumn = saveColumnAddress_;
3499  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
3500  CoinSimplexInt *  COIN_RESTRICT indexRowU = indexRowUAddress_;
3501
3502  if ( indexRowU[startColumn] == pivotRow ) {
3503    pivotElement = elementU[startColumn];
3504    otherMultiplier = elementU[startColumn + 1];
3505    otherRow = indexRowU[startColumn + 1];
3506  } else {
3507    pivotElement = elementU[startColumn + 1];
3508    otherMultiplier = elementU[startColumn];
3509    otherRow = indexRowU[startColumn];
3510  }
3511  CoinSimplexInt numberSave = numberInRow[otherRow];
3512  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
3513
3514  CoinFactorizationDouble *  COIN_RESTRICT pivotRegion = pivotRegionAddress_;
3515  pivotRegion[numberGoodU_] = pivotMultiplier;
3516  numberInColumn[pivotColumn] = 0;
3517  otherMultiplier = otherMultiplier * pivotMultiplier;
3518#ifdef SMALL_PERMUTE
3519  indexRowL[l] = fromSmallToBigRow_[otherRow];
3520#else
3521  indexRowL[l] = otherRow;
3522#endif
3523  elementL[l] = otherMultiplier;
3524  //take out of row list
3525  CoinBigIndex start = startRowU[otherRow];
3526  CoinBigIndex end = start + numberSave;
3527  CoinBigIndex where = start;
3528  CoinSimplexInt *  COIN_RESTRICT indexColumnU = indexColumnUAddress_;
3529
3530  while ( indexColumnU[where] != pivotColumn ) {
3531    where++;
3532  }                             /* endwhile */
3533  assert ( where < end );
3534#if BOTH_WAYS
3535#endif
3536  end--;
3537  indexColumnU[where] = indexColumnU[end];
3538  CoinSimplexInt numberAdded = 0;
3539  CoinSimplexInt numberDeleted = 0;
3540
3541  //pack down and move to work
3542  CoinSimplexInt j;
3543  const CoinSimplexInt *  COIN_RESTRICT nextCount = this->nextCountAddress_;
3544  CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumnAddress_;
3545  for ( j = startRow; j < endRow; j++ ) {
3546    CoinSimplexInt iColumn = indexColumnU[j];
3547
3548    if ( iColumn != pivotColumn ) {
3549      CoinBigIndex startColumn = startColumnU[iColumn];
3550      CoinSimplexInt numberInColumnIn=numberInColumn[iColumn];
3551      CoinBigIndex endColumn = startColumn + numberInColumnIn;
3552      CoinSimplexInt iRow = indexRowU[startColumn];
3553      CoinFactorizationDouble value = elementU[startColumn];
3554      CoinFactorizationDouble largest;
3555      bool foundOther = false;
3556
3557      //leave room for pivot
3558      CoinBigIndex put = startColumn + 1;
3559      CoinBigIndex positionLargest = -1;
3560      CoinFactorizationDouble thisPivotValue = 0.0;
3561      CoinFactorizationDouble otherElement = 0.0;
3562      CoinFactorizationDouble nextValue = elementU[put];;
3563      CoinSimplexInt nextIRow = indexRowU[put];
3564
3565      //compress column and find largest not updated
3566      if ( iRow != pivotRow ) {
3567        if ( iRow != otherRow ) {
3568          largest = fabs ( value );
3569          elementU[put] = value;
3570          indexRowU[put] = iRow;
3571          positionLargest = put;
3572#if BOTH_WAYS
3573#endif
3574          put++;
3575          CoinBigIndex i;
3576          for ( i = startColumn + 1; i < endColumn; i++ ) {
3577            iRow = nextIRow;
3578            value = nextValue;
3579            nextIRow = indexRowU[i + 1];
3580            nextValue = elementU[i + 1];
3581            if ( iRow != pivotRow ) {
3582              if ( iRow != otherRow ) {
3583                //keep
3584                indexRowU[put] = iRow;
3585                elementU[put] = value;;
3586#if BOTH_WAYS
3587#endif
3588                put++;
3589              } else {
3590                otherElement = value;
3591                foundOther = true;
3592              }
3593            } else {
3594              thisPivotValue = value;
3595            }
3596          }
3597        } else {
3598          otherElement = value;
3599          foundOther = true;
3600          //need to find largest
3601          largest = 0.0;
3602          CoinBigIndex i;
3603          for ( i = startColumn + 1; i < endColumn; i++ ) {
3604            iRow = nextIRow;
3605            value = nextValue;
3606            nextIRow = indexRowU[i + 1];
3607            nextValue = elementU[i + 1];
3608            if ( iRow != pivotRow ) {
3609              //keep
3610              indexRowU[put] = iRow;
3611              elementU[put] = value;;
3612#if BOTH_WAYS
3613#endif
3614              CoinFactorizationDouble absValue = fabs ( value );
3615
3616              if ( absValue > largest ) {
3617                largest = absValue;
3618                positionLargest = put;
3619              }
3620              put++;
3621            } else {
3622              thisPivotValue = value;
3623            }
3624          }
3625        }
3626      } else {
3627        //need to find largest
3628        largest = 0.0;
3629        thisPivotValue = value;
3630        CoinBigIndex i;
3631        for ( i = startColumn + 1; i < endColumn; i++ ) {
3632          iRow = nextIRow;
3633          value = nextValue;
3634          nextIRow = indexRowU[i + 1];
3635          nextValue = elementU[i + 1];
3636          if ( iRow != otherRow ) {
3637            //keep
3638            indexRowU[put] = iRow;
3639            elementU[put] = value;;
3640#if BOTH_WAYS
3641#endif
3642            CoinFactorizationDouble absValue = fabs ( value );
3643
3644            if ( absValue > largest ) {
3645              largest = absValue;
3646              positionLargest = put;
3647            }
3648            put++;
3649          } else {
3650            otherElement = value;
3651            foundOther = true;
3652          }
3653        }
3654      }
3655      //slot in pivot
3656      elementU[startColumn] = thisPivotValue;
3657#ifdef SMALL_PERMUTE
3658      indexRowU[startColumn] = realPivotRow;
3659#else
3660      indexRowU[startColumn] = pivotRow;
3661#endif
3662      //clean up counts
3663      startColumn++;
3664      numberInColumn[iColumn] = put - startColumn;
3665      numberInColumnPlus[iColumn]++;
3666      startColumnU[iColumn]++;
3667      otherElement = otherElement - thisPivotValue * otherMultiplier;
3668      CoinFactorizationDouble absValue = fabs ( otherElement );
3669
3670      if ( !TEST_LESS_THAN_TOLERANCE_REGISTER(absValue) ) {
3671        if ( !foundOther ) {
3672          //have we space
3673          saveColumn[numberAdded++] = iColumn;
3674          CoinSimplexInt next = nextColumn[iColumn];
3675          CoinBigIndex space;
3676
3677          space = startColumnU[next] - put - numberInColumnPlus[next];
3678          if ( space <= 0 ) {
3679            //getColumnSpace also moves fixed part
3680            CoinSimplexInt number = numberInColumn[iColumn];
3681
3682            if ( !getColumnSpace ( iColumn, number + 1 ) ) {
3683              return false;
3684            }
3685            //redo starts
3686            positionLargest =
3687              positionLargest + startColumnU[iColumn] - startColumn;
3688            startColumn = startColumnU[iColumn];
3689            put = startColumn + number;
3690          }
3691        }
3692        elementU[put] = otherElement;
3693        indexRowU[put] = otherRow;
3694#if BOTH_WAYS
3695#endif
3696        if ( absValue > largest ) {
3697          largest = absValue;
3698          positionLargest = put;
3699        }
3700        put++;
3701      } else {
3702        if ( foundOther ) {
3703          numberDeleted++;
3704          //take out of row list
3705          CoinBigIndex where = start;
3706
3707          while ( indexColumnU[where] != iColumn ) {
3708            where++;
3709          }                     /* endwhile */
3710          assert ( where < end );
3711#if BOTH_WAYS
3712#endif
3713          end--;
3714          indexColumnU[where] = indexColumnU[end];
3715        }
3716      }
3717      CoinSimplexInt numberInColumnOut = put - startColumn;
3718      numberInColumn[iColumn] = numberInColumnOut;
3719      //move largest
3720      if ( positionLargest >= 0 ) {
3721        value = elementU[positionLargest];
3722        iRow = indexRowU[positionLargest];
3723        elementU[positionLargest] = elementU[startColumn];
3724        indexRowU[positionLargest] = indexRowU[startColumn];
3725        elementU[startColumn] = value;
3726        indexRowU[startColumn] = iRow;
3727#if BOTH_WAYS
3728#endif
3729      }
3730      //linked list for column
3731      if ( numberInColumnIn!=numberInColumnOut&&nextCount[iColumn + numberRows_] != -2 ) {
3732        //modify linked list
3733        modifyLink ( iColumn + numberRows_, numberInColumnOut );
3734      }
3735    }
3736  }
3737  //get space for row list
3738  next = nextRow[otherRow];
3739  CoinBigIndex space;
3740
3741  space = startRowU[next] - end;
3742  totalElements_ += numberAdded - numberDeleted;
3743  CoinSimplexInt number = numberAdded + ( end - start );
3744
3745  if ( space < numberAdded ) {
3746    numberInRow[otherRow] = end - start;
3747    if ( !getRowSpace ( otherRow, number ) ) {
3748      return false;
3749    }
3750    end = startRowU[otherRow] + end - start;
3751  }
3752  // do linked lists and update counts
3753  numberInRow[otherRow] = number;
3754  if ( number != numberSave ) {
3755    modifyLink ( otherRow, number );
3756  }
3757  for ( j = 0; j < numberAdded; j++ ) {
3758#if BOTH_WAYS
3759#endif
3760    indexColumnU[end++] = saveColumn[j];
3761  }
3762  lastEntryByRowU_ = CoinMax(end,lastEntryByRowU_);
3763  assert (lastEntryByRowU_<=lengthAreaU_);
3764  //modify linked list for pivots
3765  deleteLink ( pivotRow );
3766  deleteLink ( pivotColumn + numberRows_ );
3767  return true;
3768}
3769// Delete all stuff
3770void 
3771CoinAbcTypeFactorization::almostDestructor()
3772{
3773  gutsOfDestructor(1);
3774}
3775// PreProcesses column ordered copy of basis
3776void 
3777CoinAbcTypeFactorization::preProcess ( )
3778{
3779  CoinBigIndex numberElements = startColumnUAddress_[numberRows_-1]
3780    + numberInColumnAddress_[numberRows_-1];
3781  setNumberElementsU(numberElements);
3782}
3783// Return largest
3784double 
3785CoinAbcTypeFactorization::preProcess3 ( )
3786{
3787  CoinSimplexInt * COIN_RESTRICT numberInRow = numberInRowAddress_;
3788  CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
3789  CoinSimplexInt * COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
3790  CoinBigIndex * COIN_RESTRICT startRow = startRowUAddress_;
3791  CoinBigIndex * COIN_RESTRICT startColumn = startColumnUAddress_;
3792
3793  totalElements_ = lengthU_;
3794  //links and initialize pivots
3795  CoinZeroN ( numberInColumn+numberRows_, maximumRowsExtra_ + 1 - numberRows_);
3796  CoinSimplexInt * COIN_RESTRICT lastRow = lastRowAddress_;
3797  CoinSimplexInt * COIN_RESTRICT nextRow = nextRowAddress_;
3798 
3799  CoinFillN ( pivotColumnAddress_, numberRows_, -1 );
3800  CoinZeroN ( numberInColumnPlus, maximumRowsExtra_ + 1 );
3801  CoinZeroN ( lastRow, numberRows_ );
3802  startColumn[maximumRowsExtra_] = lengthU_;
3803  lastEntryByColumnU_=lengthU_;
3804  lastEntryByRowU_=lengthU_;
3805 
3806  //do slacks first
3807  //CoinBigIndex * startColumnU = startColumnUAddress_;
3808  CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
3809  CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
3810  //CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
3811  CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
3812  CoinFactorizationDouble *  COIN_RESTRICT pivotRegion = pivotRegionAddress_;
3813  CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_;
3814  //CoinFillN(pivotRegion,numberSlacks_,static_cast<double>(SLACK_VALUE));
3815  // **************** TEMP *********************
3816  CoinFillN(pivotRegion,numberSlacks_,static_cast<CoinFactorizationDouble>(1.0));
3817  CoinZeroN(startColumnL,numberSlacks_+1);
3818  CoinZeroN(numberInColumn,numberSlacks_);
3819  CoinZeroN(numberInColumnPlus,numberSlacks_);
3820#ifdef ABC_USE_FUNCTION_POINTERS
3821  // Use row copy
3822  elementUColumnPlusAddress_ = elementRowU_.array();
3823  scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
3824#if ABC_USE_FUNCTION_POINTERS
3825  extern scatterUpdate AbcScatterLowSubtract[9];
3826#endif
3827  lastEntryByColumnUPlus_=1;
3828  //for (int i=0;i<numberRows_;i++) {
3829  //scatter[iRow].functionPointer=AbcScatterLow[0];
3830  //scatter[iRow].offset=1;
3831  //scatter[iRow].number=0;
3832  //}
3833  //extern scatterUpdate AbcScatterHigh[4];
3834#endif
3835  CoinSimplexInt * COIN_RESTRICT permute = permuteAddress_;
3836  CoinFillN(permute,numberRows_,-1);
3837  for (CoinSimplexInt iSlack = 0; iSlack < numberSlacks_;iSlack++ ) {
3838    // slack
3839    CoinSimplexInt iRow = indexRow[iSlack];
3840    lastRow[iRow] =-2;
3841    //nextRow[iRow] = numberGoodU_;     //use for permute
3842    permute[iRow]=iSlack; 
3843    numberInRow[iRow]=0;
3844    pivotColumn[iSlack] = iSlack;
3845#ifdef ABC_USE_FUNCTION_POINTERS
3846#if ABC_USE_FUNCTION_POINTERS
3847    scatter[iRow].functionPointer=AbcScatterLowSubtract[0];
3848#endif
3849    scatter[iRow].offset=0;
3850    scatter[iRow].number=0;
3851#endif
3852  }
3853  totalElements_ -= numberSlacks_;
3854  numberGoodU_=numberSlacks_;
3855  numberGoodL_=numberSlacks_;
3856#if CONVERTROW>1
3857  CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
3858  CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
3859#endif
3860#ifdef SMALL_PERMUTE
3861  CoinSimplexInt * COIN_RESTRICT fromBigToSmallRow=reinterpret_cast<CoinSimplexInt *>(saveColumn_.array());
3862  CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
3863  //move largest in column to beginning
3864  // If we always compress then can combine
3865  CoinSimplexInt *  COIN_RESTRICT tempIndex = firstCountAddress_; //saveColumnAddress_;
3866  CoinFactorizationDouble *  COIN_RESTRICT workArea = workAreaAddress_;
3867  // Compress
3868  int nSmall=0;
3869  // fromBig is only used at setup time
3870  for (int i=0;i<numberRows_;i++) {
3871    if (permute[i]>=0) {
3872      // out
3873      fromBigToSmallRow[i]=-1;
3874    } else {
3875      // in
3876      fromBigToSmallRow[i]=nSmall;
3877      numberInRow[nSmall]=numberInRow[i];
3878      lastRow[nSmall]=lastRow[i];
3879      fromSmallToBigRow_[nSmall++]=i;
3880    }
3881  }
3882  nSmall=0;
3883  for (CoinSimplexInt iColumn = numberSlacks_; iColumn < numberRows_; iColumn++ ) {
3884    CoinSimplexInt number = numberInColumn[iColumn];
3885    // use workArea and tempIndex for remaining elements
3886    CoinBigIndex first = startColumn[iColumn];
3887    int saveFirst=first;
3888    fromSmallToBigColumn_[nSmall]=iColumn;
3889    CoinBigIndex largest = -1;
3890    CoinFactorizationDouble valueLargest = -1.0;
3891    CoinSimplexInt nOther=0;
3892    CoinBigIndex end = first+number;
3893    for (CoinBigIndex k = first ; k < end; k++ ) {
3894      CoinSimplexInt iRow = indexRow[k];
3895      assert (iRow<numberRowsSmall_);
3896      CoinFactorizationDouble value = element[k];
3897      if (permute[iRow]<0) {
3898        CoinFactorizationDouble valueAbs = fabs ( value );
3899        if ( valueAbs > valueLargest ) {
3900          valueLargest = valueAbs;
3901          largest = nOther;
3902        }
3903        iRow=fromBigToSmallRow[iRow];
3904        assert (iRow>=0);
3905        tempIndex[nOther]=iRow;
3906        workArea[nOther++]=value;
3907      } else {
3908        indexRow[first] = iRow;
3909        element[first++] = value;
3910      }
3911    }
3912    CoinSimplexInt nMoved = first-saveFirst;
3913    numberInColumnPlus[nSmall]=nMoved;
3914    totalElements_ -= nMoved;
3915    startColumn[nSmall]=first;
3916    //largest
3917    if (largest>=0) {
3918      indexRow[first] = tempIndex[largest];
3919      element[first++] = workArea[largest];
3920    }
3921    for (CoinBigIndex k=0;k<nOther;k++) {
3922      if (k!=largest) {
3923        indexRow[first] = tempIndex[k];
3924        element[first++] = workArea[k];
3925      }
3926    }
3927    numberInColumn[nSmall]=first-startColumn[nSmall];
3928    nSmall++;
3929  }
3930  numberRowsSmall_=nSmall;
3931  CoinAbcMemset0(numberInRow+nSmall,numberRows_-nSmall);
3932  CoinAbcMemset0(numberInColumn+nSmall,numberRows_-nSmall);
3933#endif
3934  //row part
3935  CoinBigIndex i = 0;
3936  for (CoinSimplexInt iRow = 0; iRow < numberRowsSmall_; iRow++ ) {
3937    startRow[iRow] = i;
3938    CoinSimplexInt n=numberInRow[iRow];
3939#ifndef SMALL_PERMUTE
3940    numberInRow[iRow]=0;
3941#endif
3942    i += n;
3943  }
3944  startRow[numberRows_]=lengthAreaU_;//i;
3945  double overallLargest=1.0;
3946#ifndef SMALL_PERMUTE
3947  CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
3948  //move largest in column to beginning
3949  CoinSimplexInt *  COIN_RESTRICT tempIndex = firstCountAddress_; //saveColumnAddress_;
3950  CoinFactorizationDouble *  COIN_RESTRICT workArea = workAreaAddress_;
3951  for (CoinSimplexInt iColumn = numberSlacks_; iColumn < numberRows_; iColumn++ ) {
3952    CoinSimplexInt number = numberInColumn[iColumn];
3953    // use workArea and tempIndex for remaining elements
3954    CoinBigIndex first = startColumn[iColumn];
3955    CoinBigIndex largest = -1;
3956    CoinFactorizationDouble valueLargest = -1.0;
3957    CoinSimplexInt nOther=0;
3958    CoinBigIndex end = first+number;
3959    for (CoinBigIndex k = first ; k < end; k++ ) {
3960      CoinSimplexInt iRow = indexRow[k];
3961      assert (iRow<numberRowsSmall_);
3962      CoinFactorizationDouble value = element[k];
3963      if (!lastRow[iRow]) {
3964        numberInRow[iRow]++;
3965        CoinFactorizationDouble valueAbs = fabs ( value );
3966        if ( valueAbs > valueLargest ) {
3967          valueLargest = valueAbs;
3968          largest = nOther;
3969        }
3970        tempIndex[nOther]=iRow;
3971        workArea[nOther++]=value;
3972#if CONVERTROW<2
3973        CoinSimplexInt iLook = startRow[iRow];
3974        startRow[iRow] = iLook + 1;
3975        indexColumn[iLook] = iColumn;
3976#endif
3977#if BOTH_WAYS
3978#endif
3979      } else {
3980        indexRow[first] = iRow;
3981#if BOTH_WAYS
3982#endif
3983        element[first++] = value;
3984      }
3985    }
3986    overallLargest=CoinMax(overallLargest,valueLargest);
3987    CoinSimplexInt nMoved = first-startColumn[iColumn];
3988    numberInColumnPlus[iColumn]=nMoved;
3989    totalElements_ -= nMoved;
3990    startColumn[iColumn]=first;
3991    //largest
3992    if (largest>=0) {
3993#if CONVERTROW>1
3994      int iRow=tempIndex[largest];
3995      indexRow[first] = iRow;
3996      CoinSimplexInt iLook = startRow[iRow];
3997      startRow[iRow] = iLook + 1;
3998      indexColumn[iLook] = iColumn;
3999      convertRowToColumn[iLook]=first;
4000      convertColumnToRow[first]=iLook;
4001#else
4002      indexRow[first] = tempIndex[largest];
4003#endif
4004#if BOTH_WAYS
4005#endif
4006      element[first++] = workArea[largest];
4007    }
4008    for (CoinBigIndex k=0;k<nOther;k++) {
4009      if (k!=largest) {
4010#if CONVERTROW>1
4011      int iRow=tempIndex[k];
4012      indexRow[first] = iRow;
4013      CoinSimplexInt iLook = startRow[iRow];
4014      startRow[iRow] = iLook + 1;
4015      indexColumn[iLook] = iColumn;
4016      convertRowToColumn[iLook]=first;
4017      convertColumnToRow[first]=iLook;
4018#else
4019      indexRow[first] = tempIndex[k];
4020#endif
4021#if BOTH_WAYS
4022#endif
4023        element[first++] = workArea[k];
4024      }
4025    }
4026    numberInColumn[iColumn]=first-startColumn[iColumn];
4027  }
4028#else
4029#if CONVERTROW>1
4030  for (CoinSimplexInt iColumn = 0; iColumn < numberRowsSmall_; iColumn++ ) {
4031    CoinSimplexInt number = numberInColumn[iColumn];
4032    CoinBigIndex first = startColumn[iColumn];
4033    CoinBigIndex end = first+number;
4034    for (CoinBigIndex k = first ; k < end; k++ ) {
4035      CoinSimplexInt iRow = indexRow[k];
4036      CoinSimplexInt iLook = startRow[iRow];
4037      startRow[iRow] = iLook + 1;
4038      indexColumn[iLook] = iColumn;
4039      convertRowToColumn[iLook]=k;
4040      convertColumnToRow[k]=iLook;
4041    }
4042  }
4043#endif
4044#endif
4045  for (CoinSimplexInt iRow=numberRowsSmall_-1;iRow>0;iRow--) 
4046    startRow[iRow]=startRow[iRow-1];
4047  startRow[0]=0;
4048  // Do links (do backwards to get rows first?)
4049  CoinSimplexInt * firstCount=tempIndex;
4050  for (int i=0;i<numberRows_+2;i++) {
4051    firstCount[i]=i-numberRows_-2;
4052    //assert (&nextCount_[firstCount[i]]==firstCount+i);
4053  }
4054  startColumnLAddress_[0] = 0;  //for luck
4055  CoinSimplexInt *lastColumn = lastColumnAddress_;
4056  CoinSimplexInt *nextColumn = nextColumnAddress_;
4057#ifdef SMALL_PERMUTE
4058  // ? do I need to do slacks (and if so can I do faster)
4059  for (CoinSimplexInt iColumn = 0; iColumn <numberRowsSmall_; iColumn++ ) {
4060    lastColumn[iColumn] = iColumn - 1;
4061    nextColumn[iColumn] = iColumn + 1;
4062    CoinSimplexInt number = numberInColumn[iColumn];
4063    addLink ( iColumn + numberRows_, number );
4064  }
4065  lastColumn[maximumRowsExtra_] = numberRowsSmall_ - 1;
4066  nextColumn[maximumRowsExtra_] = 0;
4067  lastColumn[0] = maximumRowsExtra_;
4068  if (numberRowsSmall_)
4069    nextColumn[numberRowsSmall_ - 1] = maximumRowsExtra_;
4070#else
4071  lastColumn[maximumRowsExtra_] = numberRows_ - 1;
4072  for (CoinSimplexInt iColumn = 0; iColumn <numberRows_; iColumn++ ) {
4073    lastColumn[iColumn] = iColumn - 1;
4074    nextColumn[iColumn] = iColumn + 1;
4075    CoinSimplexInt number = numberInColumn[iColumn];
4076    addLink ( iColumn + numberRows_, number );
4077  }
4078  nextColumn[maximumRowsExtra_] = 0;
4079  lastColumn[0] = maximumRowsExtra_;
4080  if (numberRows_)
4081    nextColumn[numberRows_ - 1] = maximumRowsExtra_;
4082#endif
4083  startColumn[maximumRowsExtra_] = lengthU_;
4084  lastEntryByColumnU_=lengthU_;
4085  CoinSimplexInt jRow=numberRows_;
4086  for (CoinSimplexInt iRow = 0; iRow < numberRowsSmall_; iRow++ ) {
4087    CoinSimplexInt number = numberInRow[iRow];
4088    if (!lastRow[iRow]) {
4089      lastRow[iRow]=jRow;
4090      nextRow[jRow]=iRow;
4091      jRow=iRow;
4092      addLink ( iRow, number );
4093    }
4094  }
4095  // Is this why startRow[max] must be 0
4096  lastRow[numberRows_] = jRow;
4097  nextRow[jRow] = numberRows_;
4098  //checkLinks(999);
4099  return overallLargest;
4100}
4101// Does post processing on valid factorization - putting variables on correct rows
4102void 
4103CoinAbcTypeFactorization::postProcess(const CoinSimplexInt * sequence, CoinSimplexInt * pivotVariable)
4104{
4105  const CoinSimplexInt *  COIN_RESTRICT back = firstCountAddress_+numberRows_;
4106  if (pivotVariable) {
4107#ifndef NDEBUG
4108    CoinFillN(pivotVariable, numberRows_, -1);
4109#endif
4110    // Redo pivot order
4111    for (CoinSimplexInt i = 0; i < numberRows_; i++) {
4112      CoinSimplexInt k = sequence[i];
4113      CoinSimplexInt j = back[i];
4114      assert (pivotVariable[j] == -1);
4115      pivotVariable[j] = k;
4116    }
4117  }
4118  // Set up permutation vector
4119  // these arrays start off as copies of permute
4120  // (and we could use permute_ instead of pivotColumn (not back though))
4121  CoinAbcMemcpy( pivotColumnAddress_  , permuteAddress_, numberRows_ );
4122  // See if worth going sparse and when
4123  checkSparse();
4124  // Set pivotRegion in case we go off end
4125  pivotRegionAddress_[-1]=0.0;
4126  pivotRegionAddress_[numberRows_]=0.0;
4127  const CoinSimplexInt *  COIN_RESTRICT pivotLOrder = this->pivotLOrderAddress_;
4128  CoinSimplexInt *  COIN_RESTRICT pivotLinkedForwards = this->pivotLinkedForwardsAddress_;
4129  CoinSimplexInt *  COIN_RESTRICT pivotLinkedBackwards = this->pivotLinkedBackwardsAddress_;
4130  CoinSimplexInt iThis=-1;
4131  for (int i=0;i<numberSlacks_;i++) {
4132#ifndef ABC_ORDERED_FACTORIZATION
4133    CoinSimplexInt iNext=pivotLOrder[i];
4134#else
4135    CoinSimplexInt iNext=i;
4136#endif
4137    pivotLinkedForwards[iThis]=iNext;
4138    iThis=iNext;
4139  }
4140  lastSlack_ = iThis;
4141  for (int i=numberSlacks_;i<numberRows_;i++) {
4142#ifndef ABC_ORDERED_FACTORIZATION
4143    CoinSimplexInt iNext=pivotLOrder[i];
4144#else
4145    CoinSimplexInt iNext=i;
4146#endif
4147    pivotLinkedForwards[iThis]=iNext;
4148    iThis=iNext;
4149  }
4150  pivotLinkedForwards[iThis]=numberRows_;
4151  pivotLinkedForwards[numberRows_]=-1;
4152  iThis=numberRows_;
4153  for (int i=numberRows_-1;i>=0;i--) {
4154#ifndef ABC_ORDERED_FACTORIZATION
4155    CoinSimplexInt iNext=pivotLOrder[i];
4156#else
4157    CoinSimplexInt iNext=i;
4158#endif
4159    pivotLinkedBackwards[iThis]=iNext;
4160    iThis=iNext;
4161  }
4162  pivotLinkedBackwards[iThis]=-1;
4163  pivotLinkedBackwards[-1]=-2;
4164#if 0 //ndef NDEBUG
4165  {
4166    char * test = new char[numberRows_];
4167    memset(test,0,numberRows_);
4168    CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
4169    CoinSimplexInt *  COIN_RESTRICT indexRowL = indexRowLAddress_;
4170    const CoinSimplexInt * COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
4171    //CoinSimplexInt *  COIN_RESTRICT pivotLOrder = firstCountAddress_;
4172    //const CoinSimplexInt * COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
4173    for (int i=0;i<numberRows_;i++) {
4174      //CoinSimplexInt pivotOrder=permuteAddress_[i];
4175#ifndef ABC_ORDERED_FACTORIZATION
4176      CoinSimplexInt iRow=pivotLOrder[i];
4177#else
4178      CoinSimplexInt iRow=i;
4179#endif
4180      assert(iRow>=0&&iRow<numberRows_);
4181      test[iRow]=1;
4182      CoinBigIndex start=startColumnL[i];
4183      CoinBigIndex end=startColumnL[i+1];
4184      for (CoinBigIndex j = start;j < end; j ++ ) {
4185        CoinSimplexInt jRow = indexRowL[j];
4186        assert (test[jRow]==0);
4187      }
4188    }
4189    delete [] test;
4190  }
4191#endif
4192  //memcpy(pivotColumnAddress_,pivotLOrder,numberRows_*sizeof(CoinSimplexInt));
4193  if (gotRCopy()) {
4194    //if (numberInColumnPlusAddress_) {
4195    CoinBigIndex *  COIN_RESTRICT startR = startColumnRAddress_ + maximumPivots_+1;
4196    // do I need to zero
4197    assert(startR+maximumRowsExtra_+1<=firstCountAddress_+CoinMax(5*numberRows_,4*numberRows_+2*maximumPivots_+4)+2 );
4198    CoinZeroN (startR,(maximumRowsExtra_+1));
4199  }
4200  goSparse2();
4201#ifndef ABC_USE_FUNCTION_POINTERS
4202  CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
4203  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnUAddress_;
4204  CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
4205  CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
4206#else
4207  scatterStruct * COIN_RESTRICT scatterPointer = scatterUColumn();
4208  CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
4209#endif
4210  //lastEntryByRowU_ = totalElements_;
4211#if ABC_SMALL<2
4212  CoinSimplexInt * COIN_RESTRICT indexRowU = indexRowUAddress_;
4213  if (gotUCopy()) {
4214    CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
4215    CoinBigIndex *  COIN_RESTRICT startRow = startRowUAddress_;
4216    CoinSimplexInt *  COIN_RESTRICT nextRow = nextRowAddress_;
4217    CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
4218    CoinSimplexInt * COIN_RESTRICT indexColumnU = indexColumnUAddress_;
4219#ifndef ABC_USE_FUNCTION_POINTERS
4220    for ( CoinSimplexInt i = 0; i < numberU_; i++ ) {
4221      CoinBigIndex start = startColumnU[i];
4222      CoinBigIndex end = start + numberInColumn[i];
4223      CoinBigIndex j;
4224      // only needed if row copy
4225      for ( j = start; j < end; j++ ) {
4226        CoinSimplexInt iRow = indexRowU[j];
4227        numberInRow[iRow]++;
4228      }
4229    }
4230#else
4231    for ( CoinSimplexInt i = 0; i < numberU_; i++ ) {
4232      scatterStruct & COIN_RESTRICT scatter = scatterPointer[i];
4233      CoinSimplexInt number = scatter.number;
4234      const CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter.offset;
4235      const CoinBigIndex * COIN_RESTRICT indices = reinterpret_cast<const CoinBigIndex *>(area+number);
4236      for (int j = 0; j < number; j++ ) {
4237        CoinSimplexInt iRow = indices[j];
4238        assert (iRow>=0&&iRow<numberRows_);
4239        numberInRow[iRow]++;
4240      }
4241    }
4242#endif
4243    CoinFactorizationDouble * COIN_RESTRICT elementRow = elementRowUAddress_;
4244    CoinBigIndex j = 0;
4245   
4246    CoinSimplexInt iRow;
4247    for ( iRow = 0; iRow < numberRows_; iRow++ ) {
4248      startRow[iRow] = j;
4249      j += numberInRow[iRow];
4250    }
4251    lastEntryByRowU_ = j;
4252   
4253#if CONVERTROW
4254    CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
4255#if CONVERTROW>2
4256    CoinBigIndex * COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
4257#endif
4258#endif
4259#ifndef ABC_USE_FUNCTION_POINTERS
4260    CoinZeroN ( numberInRow, numberRows_+1 );
4261    for (CoinSimplexInt i = 0; i < numberRows_; i++ ) {
4262      CoinBigIndex start = startColumnU[i];
4263      CoinBigIndex end = start + numberInColumn[i];
4264      CoinFactorizationDouble pivotValue = pivotRegion[i];
4265      CoinBigIndex j;
4266      for ( j = start; j < end; j++ ) {
4267        CoinSimplexInt iRow = indexRowU[j];
4268        CoinSimplexInt iLook = numberInRow[iRow];
4269        numberInRow[iRow] = iLook + 1;
4270        CoinBigIndex k = startRow[iRow] + iLook;
4271        indexColumnU[k] = i;
4272#if CONVERTROW
4273        convertRowToColumn[k] = j-start;
4274#if CONVERTROW>2
4275        convertColumnToRow[j] = k;
4276#endif
4277#endif
4278        //multiply by pivot
4279        elementU[j] *= pivotValue;
4280        elementRow[k] = elementU[j];
4281      }
4282    }
4283#else
4284    for (CoinSimplexInt i = 0; i < numberRows_; i++ ) {
4285      scatterStruct & COIN_RESTRICT scatter = scatterPointer[i];
4286      CoinSimplexInt number = scatter.number;
4287      const CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+scatter.offset;
4288      const CoinBigIndex * COIN_RESTRICT indices = reinterpret_cast<const CoinBigIndex *>(area+number);
4289      for (int j = 0; j < number; j++ ) {
4290        CoinSimplexInt iRow = indices[j];
4291        CoinBigIndex k = startRow[iRow];
4292        startRow[iRow]=k+1;
4293        indexColumnU[k] = i;
4294#if CONVERTROW
4295        convertRowToColumn[k] = j;
4296#if CONVERTROW>2
4297        convertColumnToRow[j+startColumnU[i]] = k;
4298        abort(); // fix
4299#endif
4300#endif
4301        //CoinFactorizationDouble pivotValue = pivotRegion[i];
4302        //multiply by pivot
4303        //elementU[j+startColumnU[i]] *= pivotValue;
4304        elementRow[k] = area[j];
4305      }
4306    }
4307#endif
4308    //CoinSimplexInt *  COIN_RESTRICT nextRow = nextRowAddress_;
4309    //CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
4310    for (int j = numberRows_-1; j >=0; j-- ) {
4311      lastRow[j] = j - 1;
4312      nextRow[j] = j + 1;
4313#ifdef ABC_USE_FUNCTION_POINTERS
4314      startRow[j+1]=startRow[j];
4315#endif
4316    }
4317    startRow[0]=0;
4318    nextRow[numberRows_ - 1] = numberRows_;
4319    lastRow[numberRows_] = numberRows_ - 1;
4320    nextRow[numberRows_] = 0;
4321    lastRow[0] = numberRows_;
4322    //memcpy(firstCountAddress_+numberRows_+1,pivotLinkedBackwards-1,(numberRows_+2)*sizeof(int));
4323  } else {
4324#endif
4325#ifndef ABC_USE_FUNCTION_POINTERS
4326    for (CoinSimplexInt i = 0; i < numberRows_; i++ ) {
4327      CoinBigIndex start = startColumnU[i];
4328      CoinBigIndex end = start + numberInColumn[i];
4329      CoinFactorizationDouble pivotValue = pivotRegion[i];
4330      CoinBigIndex j;
4331      for ( j = start; j < end; j++ ) {
4332        //multiply by pivot
4333        elementU[j] *= pivotValue;
4334      }
4335    }
4336#endif
4337#if ABC_SMALL<2
4338  }
4339#endif
4340#ifdef ABC_USE_FUNCTION_POINTERS
4341#if 0
4342  {
4343    CoinFactorizationDouble *  COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
4344    CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
4345    CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
4346    CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
4347    CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
4348    scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
4349    extern scatterUpdate AbcScatterLowSubtract[9];
4350    extern scatterUpdate AbcScatterHighSubtract[4];
4351    for (int iRow=0;iRow<numberRows_;iRow++) {
4352      int number=scatter[iRow].number;
4353      CoinBigIndex offset = scatter[iRow].offset;
4354      CoinFactorizationDouble * COIN_RESTRICT area = elementUColumnPlus+offset;
4355      CoinSimplexInt * COIN_RESTRICT indices = reinterpret_cast<CoinSimplexInt *>(area+number);
4356      CoinSimplexInt * COIN_RESTRICT indices2 = indexRow+startColumnU[iRow];
4357      CoinFactorizationDouble * COIN_RESTRICT area2 = element+startColumnU[iRow];
4358      assert (number==numberInColumn[iRow]);
4359      if (number<9)
4360        assert (scatter[iRow].functionPointer==AbcScatterLowSubtract[number]);
4361      else
4362        assert (scatter[iRow].functionPointer==AbcScatterHighSubtract[number&3]);
4363      for (int i=0;i<number;i++) {
4364        assert (indices[i]==indices2[i]);
4365        assert (area[i]==area2[i]);
4366      }
4367    }
4368  }
4369#endif
4370#endif
4371}
4372// Makes a non-singular basis by replacing variables
4373void 
4374CoinAbcTypeFactorization::makeNonSingular(CoinSimplexInt *  COIN_RESTRICT sequence)
4375{
4376  // Replace bad ones by correct slack
4377  CoinSimplexInt *  COIN_RESTRICT workArea = indexRowUAddress_;
4378  for (CoinSimplexInt i=0;i<numberRows_;i++) 
4379    workArea[i]=-1;
4380  const CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_;
4381  const CoinSimplexInt *  COIN_RESTRICT permute = permuteAddress_;
4382  for ( CoinSimplexInt i=0;i<numberGoodU_;i++) {
4383    CoinSimplexInt iOriginal = pivotColumn[i];
4384    assert (iOriginal>=0);
4385    workArea[iOriginal]=i;
4386  }
4387  CoinSimplexInt iRow=0;
4388  for (CoinSimplexInt i=0;i<numberRows_;i++) {
4389    if (workArea[i]==-1) {
4390      while (permute[iRow]>=0)
4391        iRow++;
4392      assert (iRow<numberRows_);
4393      // Put slack in basis
4394      sequence[i]=iRow;
4395      iRow++;
4396    }
4397  }
4398}
4399
4400//  show_self.  Debug show object
4401void
4402CoinAbcTypeFactorization::show_self (  ) const
4403{
4404  CoinSimplexInt i;
4405
4406  const CoinSimplexInt * pivotColumn = pivotColumnAddress_;
4407  for ( i = 0; i < numberRows_; i++ ) {
4408    std::cout << "r " << i << " " << pivotColumn[i];
4409    std::cout<< " " << permuteAddress_[i];
4410    std::cout<< " " << pivotRegionAddress_[i];
4411    std::cout << std::endl;
4412  }
4413  for ( i = 0; i < numberRows_; i++ ) {
4414    std::cout << "u " << i << " " << numberInColumnAddress_[i] << std::endl;
4415    CoinSimplexInt j;
4416    CoinSort_2(indexRowUAddress_+startColumnUAddress_[i],
4417               indexRowUAddress_+startColumnUAddress_[i]+numberInColumnAddress_[i],
4418               elementUAddress_+startColumnUAddress_[i]);
4419    for ( j = startColumnUAddress_[i]; j < startColumnUAddress_[i] + numberInColumnAddress_[i];
4420          j++ ) {
4421      assert (indexRowUAddress_[j]>=0&&indexRowUAddress_[j]<numberRows_);
4422      assert (elementUAddress_[j]>-1.0e50&&elementUAddress_[j]<1.0e50);
4423      std::cout << indexRowUAddress_[j] << " " << elementUAddress_[j] << std::endl;
4424    }
4425  }
4426  for ( i = 0; i < numberRows_; i++ ) {
4427    std::cout << "l " << i << " " << startColumnLAddress_[i + 1] -
4428      startColumnLAddress_[i] << std::endl;
4429    CoinSort_2(indexRowLAddress_+startColumnLAddress_[i],
4430               indexRowLAddress_+startColumnLAddress_[i+1],
4431               elementLAddress_+startColumnLAddress_[i]);
4432    CoinSimplexInt j;
4433    for ( j = startColumnLAddress_[i]; j < startColumnLAddress_[i + 1]; j++ ) {
4434      std::cout << indexRowLAddress_[j] << " " << elementLAddress_[j] << std::endl;
4435    }
4436  }
4437
4438}
4439//  sort so can compare
4440void
4441CoinAbcTypeFactorization::sort (  ) const
4442{
4443  CoinSimplexInt i;
4444
4445  for ( i = 0; i < numberRows_; i++ ) {
4446    CoinSort_2(indexRowUAddress_+startColumnUAddress_[i],
4447               indexRowUAddress_+startColumnUAddress_[i]+numberInColumnAddress_[i],
4448               elementUAddress_+startColumnUAddress_[i]);
4449  }
4450  for ( i = 0; i < numberRows_; i++ ) {
4451    CoinSort_2(indexRowLAddress_+startColumnLAddress_[i],
4452               indexRowLAddress_+startColumnLAddress_[i+1],
4453               elementLAddress_+startColumnLAddress_[i]);
4454  }
4455
4456}
4457#ifdef CHECK_LINKS
4458//static int ixxxxxx=0;
4459void
4460CoinAbcTypeFactorization::checkLinks(int type)
4461{
4462  //if (type<1)
4463  return;
4464#if 0 //ABC_SMALL<2
4465  if (numberRowsSmall_!=5698)
4466    return;
4467  ixxxxxx++;
4468  if (fromSmallToBigRow_[3195]==3581) {
4469    if (permuteAddress_[3581]>=0) {
4470      assert (numberInRowAddress_[3195]==0);
4471      CoinSimplexInt * numberInColumn = numberInColumnAddress_;
4472      CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
4473      CoinBigIndex *  COIN_RESTRICT startColumn = startColumnUAddress_;
4474      for (int i=startColumn[1523];i<startColumn[1523]+numberInColumn[1523];i++)
4475        assert (indexRow[i]!=3195);
4476    }
4477    CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
4478    CoinSimplexInt *  COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
4479    CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumnAddress_;
4480    CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
4481    CoinSimplexInt iColumn = nextColumn[maximumRowsExtra_];
4482    //CoinBigIndex put = startColumnU[iColumn] - numberInColumnPlus[iColumn];
4483    //CoinBigIndex putEnd = startColumnU[iColumn] + numberInColumn[iColumn];
4484    CoinBigIndex putEnd = -1;
4485    while ( iColumn != maximumRowsExtra_ ) {
4486      CoinBigIndex get;
4487      CoinBigIndex getEnd;
4488      get = startColumnU[iColumn] - numberInColumnPlus[iColumn];
4489      getEnd = startColumnU[iColumn] + numberInColumn[iColumn];
4490      assert (putEnd<=get);
4491      putEnd=getEnd;
4492      iColumn = nextColumn[iColumn];
4493    }
4494  }
4495  if (ixxxxxx>=7615900) {
4496    printf("trouble ahead\n");
4497  }
4498  return;
4499  if (ixxxxxx<406245)
4500    return;
4501#endif
4502  CoinSimplexInt * COIN_RESTRICT nextCount = this->nextCount();
4503  CoinSimplexInt * COIN_RESTRICT firstCount = this->firstCount();
4504  CoinSimplexInt * COIN_RESTRICT lastCount = this->lastCount();
4505  CoinSimplexInt * numberInColumn = numberInColumnAddress_;
4506  CoinSimplexInt * numberInRow = numberInRowAddress_;
4507  if (type) {
4508    int nBad=0;
4509    CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
4510    CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
4511    CoinBigIndex *  COIN_RESTRICT startColumn = startColumnUAddress_;
4512    CoinBigIndex * COIN_RESTRICT startRow = startRowUAddress_;
4513    for (int i=0;i<numberRowsSmall_;i++) {
4514      assert (startRow[i]+numberInRow[i]<=lastEntryByRowU_);
4515    }
4516    return;
4517    CoinSimplexInt * COIN_RESTRICT nextRow = nextRow_.array();
4518    CoinSimplexInt * COIN_RESTRICT lastRow = lastRow_.array() ;
4519    int * counts = new int [3*numberRows_+2];
4520    int * list1=counts+numberRows_;
4521    int * list2=list1+numberRows_+1;
4522    {
4523      int lastStart=-1;
4524      int iRow=nextRow[numberRows_];
4525      memset(list1,0,numberRowsSmall_*sizeof(int));
4526      while (iRow>=0&&iRow<numberRowsSmall_) {
4527        int iStart=startRow[iRow];
4528        assert (iStart>=lastStart);
4529        lastStart=iStart;
4530        list1[iRow]=1;
4531        iRow=nextRow[iRow];
4532      }
4533      for (int i=0;i<numberRowsSmall_;i++) { 
4534        assert (list1[i]||!numberInRow[i]);
4535        assert (!list1[i]||numberInRow[i]);
4536      }
4537      memset(list1,0,numberRowsSmall_*sizeof(int));
4538#if 1
4539      CoinSimplexInt * COIN_RESTRICT nextColumn = nextColumn_.array();
4540      CoinSimplexInt * COIN_RESTRICT lastColumn = lastColumn_.array() ;
4541      lastStart=-1;
4542      int iColumn=nextColumn[maximumRowsExtra_];
4543      int nDone=0;
4544      while (iColumn>=0&&iColumn<numberRowsSmall_) {
4545        nDone++;
4546        int iStart=startColumn[iColumn];
4547        assert (iStart>=lastStart);
4548        lastStart=iStart;
4549        list1[iColumn]=1;
4550        iColumn=nextColumn[iColumn];
4551      }
4552      int nEnd=0;
4553      for (int i=0;i<numberRowsSmall_;i++) {
4554        if (numberInColumn[i]&&nEnd==maximumRowsExtra_)
4555          nEnd++;
4556      }
4557      if (nEnd>1) {
4558        for (int i=0;i<numberRowsSmall_;i++) {
4559          if (numberInColumn[i]&&nEnd==maximumRowsExtra_)
4560            printf("%d points to end\n",i);
4561        }
4562        assert (nEnd<=1);
4563      }
4564      //printf("rows %d goodU %d done %d iColumn %d\n",numberRowsSmall_,numberGoodU_,nDone,iColumn);
4565#endif
4566#if 0
4567      for (int i=0;i<numberRowsSmall_;i++) {
4568        assert (list1[i]||!numberInColumn[i]);
4569        assert (!list1[i]||numberInColumn[i]);
4570        //assert (list1[i]||!numberInColumn[fromSmallToBigColumn[i]]);
4571        //assert (!list1[i]||numberInColumn[fromSmallToBigColumn[i]]);
4572      }
4573#endif
4574    }
4575    CoinAbcMemset0(counts,numberRows_);
4576    for (int iColumn=0;iColumn<numberRows_;iColumn++) {
4577      for (CoinBigIndex j=startColumn[iColumn];j<startColumn[iColumn]+numberInColumn[iColumn];j++) {
4578        int iRow=indexRow[j];
4579        counts[iRow]++;
4580      }
4581    }
4582    int nTotal1=0;
4583    for (int i=0;i<numberRows_;i++) {
4584      if (numberInRow[i]!=counts[i]) {
4585        printf("Row %d numberInRow %d counts %d\n",i,numberInRow[i],counts[i]);
4586        int n1=0;
4587        int n2=0;
4588        for (int j=startRow[i];j<startRow[i]+numberInRow[i];j++)
4589          list1[n1++]=indexColumn[j];
4590        for (int iColumn=0;iColumn<numberRows_;iColumn++) {
4591          for (CoinBigIndex j=startColumn[iColumn];j<startColumn[iColumn]+numberInColumn[iColumn];j++) {
4592            int iRow=indexRow[j];
4593            if (i==iRow)
4594              list2[n2++]=iColumn;
4595          }
4596        }
4597        std::sort(list1,list1+n1);
4598        std::sort(list2,list2+n2);
4599        list1[n1]=numberRows_;
4600        list2[n2]=numberRows_;
4601        n1=1;
4602        n2=1;
4603        int iLook1=list1[0];
4604        int iLook2=list2[0];
4605        int xx=0;
4606        while(iLook1!=numberRows_||iLook2!=numberRows_) {
4607          if (iLook1==iLook2) {
4608            printf(" (rc %d)",iLook1);
4609            iLook1=list1[n1++];
4610            iLook2=list2[n2++];
4611          } else if (iLook1<iLook2) {
4612            printf(" (r %d)",iLook1);
4613            iLook1=list1[n1++];
4614          } else {
4615            printf(" (c %d)",iLook2);
4616            iLook2=list2[n2++];
4617          }
4618          xx++;
4619          if (xx==10) {
4620            printf("\n");
4621            xx=0;
4622          }
4623        }
4624        if (xx)
4625          printf("\n");
4626        nBad++;
4627      }
4628      nTotal1+=counts[i];
4629    }
4630    CoinAbcMemset0(counts,numberRows_);
4631    int nTotal2=0;
4632    for (int iRow=0;iRow<numberRows_;iRow++) {
4633      for (CoinBigIndex j=startRow[iRow];j<startRow[iRow]+numberInRow[iRow];j++) {
4634        int iColumn=indexColumn[j];
4635        counts[iColumn]++;
4636      }
4637    }
4638    for (int i=0;i<numberRows_;i++) {
4639      if (numberInColumn[i]!=counts[i]) {
4640        printf("Column %d numberInColumn %d counts %d\n",i,numberInColumn[i],counts[i]);
4641        int n1=0;
4642        int n2=0;
4643        for (int j=startColumn[i];j<startColumn[i]+numberInColumn[i];j++)
4644          list1[n1++]=indexRow[j];
4645        for (int iRow=0;iRow<numberRows_;iRow++) {
4646          for (CoinBigIndex j=startRow[iRow];j<startRow[iRow]+numberInRow[iRow];j++) {
4647            int iColumn=indexColumn[j];
4648            if (i==iColumn)
4649              list2[n2++]=iRow;
4650          }
4651        }
4652        std::sort(list1,list1+n1);
4653        std::sort(list2,list2+n2);
4654        list1[n1]=numberRows_;
4655        list2[n2]=numberRows_;
4656        n1=1;
4657        n2=1;
4658        int iLook1=list1[0];
4659        int iLook2=list2[0];
4660        int xx=0;
4661        while(iLook1!=numberRows_||iLook2!=numberRows_) {
4662          if (iLook1==iLook2) {
4663            printf(" (rc %d)",iLook1);
4664            iLook1=list1[n1++];
4665            iLook2=list2[n2++];
4666          } else if (iLook1<iLook2) {
4667            printf(" (c %d)",iLook1);
4668            iLook1=list1[n1++];
4669          } else {
4670            printf(" (r %d)",iLook2);
4671            iLook2=list2[n2++];
4672          }
4673          xx++;
4674          if (xx==10) {
4675            printf("\n");
4676            xx=0;
4677          }
4678        }
4679        if (xx)
4680          printf("\n");
4681        nBad++;
4682      }
4683      nTotal2+=counts[i];
4684    }
4685    delete [] counts;
4686    assert (!nBad);
4687  }
4688  char * mark = new char [2*numberRows_];
4689  memset(mark,0,2*numberRows_);
4690  for (int iCount=0;iCount<numberRowsSmall_+1;iCount++) {
4691    int next=firstCount[iCount];
4692#define VERBOSE 0
4693#if VERBOSE
4694    int count=0;
4695#endif
4696    while (next>=0) {
4697#if VERBOSE
4698      count++;
4699      if (iCount<=VERBOSE) {
4700        if (next<numberRows_)
4701#ifdef SMALL_PERMUTE
4702          printf("R%d\n",fromSmallToBigRow_[next]);
4703#else
4704          printf("R%d\n",next);
4705#endif
4706        else
4707          printf("C%d\n",next-numberRows_);
4708      }
4709#endif
4710      assert (!mark[next]);
4711      mark[next]=1;
4712      if (next<numberRows_)
4713        assert (numberInRow[next]==iCount);
4714      else
4715        assert (numberInColumn[next-numberRows_]==iCount);
4716      int next2=next;
4717      assert (next!=nextCount[next]);
4718      next=nextCount[next];
4719      if (next>=0)
4720        assert (next2==lastCount[next]);
4721    }
4722#if VERBOSE
4723    if (count)
4724      printf("%d items have count %d\n",count,iCount); 
4725#endif
4726  }
4727  delete [] mark;
4728}
4729#endif
4730#endif
Note: See TracBrowser for help on using the repository browser.