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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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