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

Last change on this file since 2271 was 2146, checked in by forrest, 4 years ago

stuff from stable and start aboca_lite

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