source: trunk/Clp/src/CoinAbcBaseFactorization3.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: 83.9 KB
Line 
1/* $Id: CoinAbcBaseFactorization3.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
13#include "CoinIndexedVector.hpp"
14#include "CoinHelperFunctions.hpp"
15#include "CoinAbcHelperFunctions.hpp"
16#if CILK_CONFLICT > 0
17// for conflicts
18extern int cilk_conflict;
19#endif
20//:class CoinAbcTypeFactorization.  Deals with Factorization and Updates
21
22/* Updates one column for dual steepest edge weights (FTRAN) */
23void CoinAbcTypeFactorization::updateWeights(CoinIndexedVector &regionSparse) const
24{
25  toLongArray(&regionSparse, 1);
26#ifdef ABC_ORDERED_FACTORIZATION
27  // Permute in for Ftran
28  permuteInForFtran(regionSparse);
29#endif
30  //  ******* L
31  updateColumnL(&regionSparse
32#if ABC_SMALL < 2
33    ,
34    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
35#endif
36#if ABC_PARALLEL
37    // can re-use 0 which would have been used for btran
38    ,
39    0
40#endif
41  );
42  //row bits here
43  updateColumnR(&regionSparse
44#if ABC_SMALL < 2
45    ,
46    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
47#endif
48#if ABC_PARALLEL
49      ,
50    0
51#endif
52  );
53  //update counts
54  //  ******* U
55  updateColumnU(&regionSparse
56#if ABC_SMALL < 2
57    ,
58    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
59#endif
60#if ABC_PARALLEL
61      ,
62    0
63#endif
64  );
65#if ABC_DEBUG
66  regionSparse.checkClean();
67#endif
68}
69CoinSimplexInt CoinAbcTypeFactorization::updateColumn(CoinIndexedVector &regionSparse)
70  const
71{
72  toLongArray(&regionSparse, 0);
73#ifdef ABC_ORDERED_FACTORIZATION
74  // Permute in for Ftran
75  permuteInForFtran(regionSparse);
76#endif
77  //  ******* L
78  updateColumnL(&regionSparse
79#if ABC_SMALL < 2
80    ,
81    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
82#endif
83  );
84  //row bits here
85  updateColumnR(&regionSparse
86#if ABC_SMALL < 2
87    ,
88    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
89#endif
90  );
91  //update counts
92  //  ******* U
93  updateColumnU(&regionSparse
94#if ABC_SMALL < 2
95    ,
96    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
97#endif
98  );
99#if ABC_DEBUG
100  regionSparse.checkClean();
101#endif
102  return regionSparse.getNumElements();
103}
104/* Updates one full column (FTRAN) */
105void CoinAbcTypeFactorization::updateFullColumn(CoinIndexedVector &regionSparse) const
106{
107#ifndef ABC_ORDERED_FACTORIZATION
108  regionSparse.setNumElements(0);
109  regionSparse.scan(0, numberRows_);
110#else
111  permuteInForFtran(regionSparse, true);
112#endif
113  if (regionSparse.getNumElements()) {
114    toLongArray(&regionSparse, 1);
115    //  ******* L
116    updateColumnL(&regionSparse
117#if ABC_SMALL < 2
118      ,
119      reinterpret_cast< CoinAbcStatistics & >(ftranFullCountInput_)
120#endif
121#if ABC_PARALLEL
122        ,
123      1
124#endif
125    );
126    //row bits here
127    updateColumnR(&regionSparse
128#if ABC_SMALL < 2
129      ,
130      reinterpret_cast< CoinAbcStatistics & >(ftranFullCountInput_)
131#endif
132#if ABC_PARALLEL
133        ,
134      1
135#endif
136    );
137    //update counts
138    //  ******* U
139    updateColumnU(&regionSparse
140#if ABC_SMALL < 2
141      ,
142      reinterpret_cast< CoinAbcStatistics & >(ftranFullCountInput_)
143#endif
144#if ABC_PARALLEL
145        ,
146      1
147#endif
148    );
149    fromLongArray(1);
150#if ABC_DEBUG
151    regionSparse.checkClean();
152#endif
153  }
154}
155// move stuff like this into CoinAbcHelperFunctions.hpp
156inline void multiplyIndexed(CoinSimplexInt number,
157  const CoinFactorizationDouble *COIN_RESTRICT multiplier,
158  const CoinSimplexInt *COIN_RESTRICT thisIndex,
159  CoinFactorizationDouble *COIN_RESTRICT region)
160{
161  for (CoinSimplexInt i = 0; i < number; i++) {
162    CoinSimplexInt iRow = thisIndex[i];
163    CoinSimplexDouble value = region[iRow];
164    value *= multiplier[iRow];
165    region[iRow] = value;
166  }
167}
168/* Updates one full column (BTRAN) */
169void CoinAbcTypeFactorization::updateFullColumnTranspose(CoinIndexedVector &regionSparse) const
170{
171  int numberNonZero = 0;
172  // Should pass in statistics
173  toLongArray(&regionSparse, 0);
174  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
175  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse.getIndices();
176#ifndef ABC_ORDERED_FACTORIZATION
177  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
178  for (int iRow = 0; iRow < numberRows_; iRow++) {
179    double value = region[iRow];
180    if (value) {
181      region[iRow] = value * pivotRegion[iRow];
182      regionIndex[numberNonZero++] = iRow;
183    }
184  }
185  regionSparse.setNumElements(numberNonZero);
186  if (!numberNonZero)
187    return;
188    //regionSparse.setNumElements(0);
189    //regionSparse.scan(0,numberRows_);
190#else
191  permuteInForBtranAndMultiply(regionSparse, true);
192  if (!regionSparse.getNumElements()) {
193    permuteOutForBtran(regionSparse);
194    return;
195  }
196#endif
197
198  //  ******* U
199  // Apply pivot region - could be combined for speed
200  // Can only combine/move inside vector for sparse
201  CoinSimplexInt smallestIndex = pivotLinkedForwardsAddress_[-1];
202#if ABC_SMALL < 2
203  // copy of code inside transposeU
204  bool goSparse = false;
205#else
206#define goSparse false
207#endif
208#if ABC_SMALL < 2
209  // Guess at number at end
210  if (gotUCopy()) {
211    assert(btranFullAverageAfterU_);
212    CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * btranFullAverageAfterU_ * twiddleBtranFullFactor1());
213    if (newNumber < sparseThreshold_)
214      goSparse = true;
215  }
216#endif
217  if (numberNonZero < 40 && (numberNonZero << 4) < numberRows_ && !goSparse) {
218    CoinSimplexInt *COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
219    CoinSimplexInt smallest = numberRowsExtra_;
220    for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
221      CoinSimplexInt iRow = regionIndex[j];
222      if (pivotRowForward[iRow] < smallest) {
223        smallest = pivotRowForward[iRow];
224        smallestIndex = iRow;
225      }
226    }
227  }
228  updateColumnTransposeU(&regionSparse, smallestIndex
229#if ABC_SMALL < 2
230    ,
231    reinterpret_cast< CoinAbcStatistics & >(btranFullCountInput_)
232#endif
233#if ABC_PARALLEL
234      ,
235    0
236#endif
237  );
238  //row bits here
239  updateColumnTransposeR(&regionSparse
240#if ABC_SMALL < 2
241    ,
242    reinterpret_cast< CoinAbcStatistics & >(btranFullCountInput_)
243#endif
244  );
245  //  ******* L
246  updateColumnTransposeL(&regionSparse
247#if ABC_SMALL < 2
248    ,
249    reinterpret_cast< CoinAbcStatistics & >(btranFullCountInput_)
250#endif
251#if ABC_PARALLEL
252      ,
253    0
254#endif
255  );
256  fromLongArray(0);
257#if ABC_SMALL < 3
258#ifdef ABC_ORDERED_FACTORIZATION
259  permuteOutForBtran(regionSparse);
260#endif
261#if ABC_DEBUG
262  regionSparse.checkClean();
263#endif
264#else
265  numberNonZero = 0;
266  for (CoinSimplexInt i = 0; i < numberRows_; i++) {
267    CoinExponent expValue = ABC_EXPONENT(region[i]);
268    if (expValue) {
269      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
270        regionIndex[numberNonZero++] = i;
271      } else {
272        region[i] = 0.0;
273      }
274    }
275  }
276  regionSparse.setNumElements(numberNonZero);
277#endif
278}
279//  updateColumnL.  Updates part of column (FTRANL)
280void CoinAbcTypeFactorization::updateColumnL(CoinIndexedVector *regionSparse
281#if ABC_SMALL < 2
282  ,
283  CoinAbcStatistics &statistics
284#endif
285#if ABC_PARALLEL
286  ,
287  int whichSparse
288#endif
289  ) const
290{
291#if CILK_CONFLICT > 0
292#if ABC_PARALLEL
293  // for conflicts
294#if CILK_CONFLICT > 1
295  printf("file %s line %d which %d\n", __FILE__, __LINE__, whichSparse);
296#endif
297  abc_assert((cilk_conflict & (1 << whichSparse)) == 0);
298  cilk_conflict |= (1 << whichSparse);
299#else
300  abc_assert((cilk_conflict & (1 << 0)) == 0);
301  cilk_conflict |= (1 << 0);
302#endif
303#endif
304#if ABC_SMALL < 2
305  CoinSimplexInt number = regionSparse->getNumElements();
306  if (factorizationStatistics()) {
307    statistics.numberCounts_++;
308    statistics.countInput_ += number;
309  }
310#endif
311  if (numberL_) {
312#if ABC_SMALL < 2
313    int goSparse;
314    // Guess at number at end
315    if (gotSparse()) {
316      double average = statistics.averageAfterL_ * twiddleFactor1S();
317      assert(average);
318      CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(number * average);
319      if (newNumber < sparseThreshold_ && (numberL_ << 2) > newNumber * 1.0 * twiddleFactor2S())
320        goSparse = 1;
321      else if (3 * newNumber < numberRows_)
322        goSparse = 0;
323      else
324        goSparse = -1;
325    } else {
326      goSparse = 0;
327    } //if(goSparse==1) goSparse=0;
328    if (!goSparse) {
329      // densish
330      updateColumnLDensish(regionSparse);
331    } else if (goSparse < 0) {
332      // densish
333      updateColumnLDense(regionSparse);
334    } else {
335      // sparse
336      updateColumnLSparse(regionSparse
337#if ABC_PARALLEL
338        ,
339        whichSparse
340#endif
341      );
342    }
343#else
344    updateColumnLDensish(regionSparse);
345#endif
346  }
347#if ABC_SMALL < 4
348#if ABC_DENSE_CODE > 0
349  if (numberDense_) {
350    instrument_do("CoinAbcFactorizationDense", 0.3 * numberDense_ * numberDense_);
351    //take off list
352    CoinSimplexInt lastSparse = numberRows_ - numberDense_;
353    CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
354    CoinFactorizationDouble *COIN_RESTRICT denseArea = denseAreaAddress_;
355    CoinFactorizationDouble *COIN_RESTRICT denseRegion = denseArea + leadingDimension_ * numberDense_;
356    CoinSimplexInt *COIN_RESTRICT densePermute = reinterpret_cast< CoinSimplexInt * >(denseRegion + FACTOR_CPU * numberDense_);
357    //for (int i=0;i<numberDense_;i++) {
358    //printf("perm %d %d\n",i,densePermute[i]);
359    //assert (densePermute[i]>=0&&densePermute[i]<numberRows_);
360    //}
361#if ABC_PARALLEL
362    if (whichSparse)
363      denseRegion += whichSparse * numberDense_;
364      //printf("PP %d %d %s region %x\n",whichSparse,__LINE__,__FILE__,denseRegion);
365#endif
366    CoinFactorizationDouble *COIN_RESTRICT densePut = denseRegion - lastSparse;
367    CoinZeroN(denseRegion, numberDense_);
368    bool doDense = false;
369#if ABC_SMALL < 3
370#if ABC_DENSE_CODE == 2
371    //short * COIN_RESTRICT forFtran = reinterpret_cast<short *>(densePermute+numberDense_)-lastSparse;
372    short *COIN_RESTRICT forFtran2 = reinterpret_cast< short * >(densePermute + numberDense_) + 2 * numberDense_;
373#else
374    const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
375#endif
376    CoinSimplexInt number = regionSparse->getNumElements();
377    CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
378    CoinSimplexInt i = 0;
379    while (i < number) {
380      CoinSimplexInt iRow = regionIndex[i];
381#if ABC_DENSE_CODE == 2
382      int kRow = forFtran2[iRow];
383      if (kRow != -1) {
384        doDense = true;
385        regionIndex[i] = regionIndex[--number];
386        denseRegion[kRow] = region[iRow];
387#else
388      CoinSimplexInt jRow = pivotLBackwardOrder[iRow];
389      if (jRow >= lastSparse) {
390        doDense = true;
391        regionIndex[i] = regionIndex[--number];
392        densePut[jRow] = region[iRow];
393#endif
394        region[iRow] = 0.0;
395      } else {
396        i++;
397      }
398    }
399#else
400    CoinSimplexInt *COIN_RESTRICT forFtran = densePermute + numberDense_ - lastSparse;
401    const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
402    for (CoinSimplexInt jRow = lastSparse; jRow < numberRows_; jRow++) {
403      CoinSimplexInt iRow = pivotLOrder[jRow];
404      if (region[iRow]) {
405        doDense = true;
406        //densePut[jRow]=region[iRow];
407#if ABC_DENSE_CODE == 2
408        int kRow = forFtran[jRow];
409        denseRegion[kRow] = region[iRow];
410#endif
411        region[iRow] = 0.0;
412      }
413    }
414#endif
415    if (doDense) {
416#if ABC_DENSE_CODE == 1
417#ifndef CLAPACK
418      char trans = 'N';
419      CoinSimplexInt ione = 1;
420      CoinSimplexInt info;
421      F77_FUNC(dgetrs, DGETRS)
422      (&trans, &numberDense_, &ione, denseArea, &leadingDimension_,
423        densePermute, denseRegion, &numberDense_, &info, 1);
424#else
425      clapack_dgetrs(CblasColMajor, CblasNoTrans, numberDense_, 1,
426        denseArea, leadingDimension_, densePermute, denseRegion, numberDense_);
427#endif
428#elif ABC_DENSE_CODE == 2
429      CoinAbcDgetrs('N', numberDense_, denseArea, denseRegion);
430#endif
431      const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
432      for (CoinSimplexInt i = lastSparse; i < numberRows_; i++) {
433        if (!TEST_LESS_THAN_TOLERANCE(densePut[i])) {
434#ifndef ABC_ORDERED_FACTORIZATION
435          CoinSimplexInt iRow = pivotLOrder[i];
436#else
437          CoinSimplexInt iRow = i;
438#endif
439          region[iRow] = densePut[i];
440#if ABC_SMALL < 3
441          regionIndex[number++] = iRow;
442#endif
443        }
444      }
445#if ABC_SMALL < 3
446      regionSparse->setNumElements(number);
447#endif
448    }
449  }
450  //printRegion(*regionSparse,"After FtranL");
451#endif
452#endif
453#if ABC_SMALL < 2
454  if (factorizationStatistics())
455    statistics.countAfterL_ += regionSparse->getNumElements();
456#endif
457#if CILK_CONFLICT > 0
458#if ABC_PARALLEL
459  // for conflicts
460  abc_assert((cilk_conflict & (1 << whichSparse)) != 0);
461  cilk_conflict &= ~(1 << whichSparse);
462#else
463  abc_assert((cilk_conflict & (1 << 0)) != 0);
464  cilk_conflict &= ~(1 << 0);
465#endif
466#endif
467}
468#define UNROLL 2
469#define INLINE_IT
470//#define INLINE_IT2
471inline void scatterUpdateInline(CoinSimplexInt number,
472  CoinFactorizationDouble pivotValue,
473  const CoinFactorizationDouble *COIN_RESTRICT thisElement,
474  const CoinSimplexInt *COIN_RESTRICT thisIndex,
475  CoinFactorizationDouble *COIN_RESTRICT region)
476{
477#if UNROLL == 0
478  for (CoinBigIndex j = number - 1; j >= 0; j--) {
479    CoinSimplexInt iRow = thisIndex[j];
480    CoinFactorizationDouble regionValue = region[iRow];
481    CoinFactorizationDouble value = thisElement[j];
482    assert(value);
483    region[iRow] = regionValue - value * pivotValue;
484  }
485#elif UNROLL == 1
486  if ((number & 1) != 0) {
487    number--;
488    CoinSimplexInt iRow = thisIndex[number];
489    CoinFactorizationDouble regionValue = region[iRow];
490    CoinFactorizationDouble value = thisElement[number];
491    region[iRow] = regionValue - value * pivotValue;
492  }
493  for (CoinBigIndex j = number - 1; j >= 0; j -= 2) {
494    CoinSimplexInt iRow0 = thisIndex[j];
495    CoinSimplexInt iRow1 = thisIndex[j - 1];
496    CoinFactorizationDouble regionValue0 = region[iRow0];
497    CoinFactorizationDouble regionValue1 = region[iRow1];
498    region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
499    region[iRow1] = regionValue1 - thisElement[j - 1] * pivotValue;
500  }
501#elif UNROLL == 2
502  CoinSimplexInt iRow0;
503  CoinSimplexInt iRow1;
504  CoinFactorizationDouble regionValue0;
505  CoinFactorizationDouble regionValue1;
506  switch (static_cast< unsigned int >(number)) {
507  case 0:
508    break;
509  case 1:
510    iRow0 = thisIndex[0];
511    regionValue0 = region[iRow0];
512    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
513    break;
514  case 2:
515    iRow0 = thisIndex[0];
516    iRow1 = thisIndex[1];
517    regionValue0 = region[iRow0];
518    regionValue1 = region[iRow1];
519    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
520    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
521    break;
522  case 3:
523    iRow0 = thisIndex[0];
524    iRow1 = thisIndex[1];
525    regionValue0 = region[iRow0];
526    regionValue1 = region[iRow1];
527    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
528    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
529    iRow0 = thisIndex[2];
530    regionValue0 = region[iRow0];
531    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
532    break;
533  case 4:
534    iRow0 = thisIndex[0];
535    iRow1 = thisIndex[1];
536    regionValue0 = region[iRow0];
537    regionValue1 = region[iRow1];
538    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
539    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
540    iRow0 = thisIndex[2];
541    iRow1 = thisIndex[3];
542    regionValue0 = region[iRow0];
543    regionValue1 = region[iRow1];
544    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
545    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
546    break;
547  case 5:
548    iRow0 = thisIndex[0];
549    iRow1 = thisIndex[1];
550    regionValue0 = region[iRow0];
551    regionValue1 = region[iRow1];
552    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
553    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
554    iRow0 = thisIndex[2];
555    iRow1 = thisIndex[3];
556    regionValue0 = region[iRow0];
557    regionValue1 = region[iRow1];
558    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
559    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
560    iRow0 = thisIndex[4];
561    regionValue0 = region[iRow0];
562    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
563    break;
564  case 6:
565    iRow0 = thisIndex[0];
566    iRow1 = thisIndex[1];
567    regionValue0 = region[iRow0];
568    regionValue1 = region[iRow1];
569    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
570    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
571    iRow0 = thisIndex[2];
572    iRow1 = thisIndex[3];
573    regionValue0 = region[iRow0];
574    regionValue1 = region[iRow1];
575    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
576    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
577    iRow0 = thisIndex[4];
578    iRow1 = thisIndex[5];
579    regionValue0 = region[iRow0];
580    regionValue1 = region[iRow1];
581    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
582    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
583    break;
584  case 7:
585    iRow0 = thisIndex[0];
586    iRow1 = thisIndex[1];
587    regionValue0 = region[iRow0];
588    regionValue1 = region[iRow1];
589    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
590    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
591    iRow0 = thisIndex[2];
592    iRow1 = thisIndex[3];
593    regionValue0 = region[iRow0];
594    regionValue1 = region[iRow1];
595    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
596    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
597    iRow0 = thisIndex[4];
598    iRow1 = thisIndex[5];
599    regionValue0 = region[iRow0];
600    regionValue1 = region[iRow1];
601    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
602    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
603    iRow0 = thisIndex[6];
604    regionValue0 = region[iRow0];
605    region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
606    break;
607  case 8:
608    iRow0 = thisIndex[0];
609    iRow1 = thisIndex[1];
610    regionValue0 = region[iRow0];
611    regionValue1 = region[iRow1];
612    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
613    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
614    iRow0 = thisIndex[2];
615    iRow1 = thisIndex[3];
616    regionValue0 = region[iRow0];
617    regionValue1 = region[iRow1];
618    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
619    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
620    iRow0 = thisIndex[4];
621    iRow1 = thisIndex[5];
622    regionValue0 = region[iRow0];
623    regionValue1 = region[iRow1];
624    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
625    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
626    iRow0 = thisIndex[6];
627    iRow1 = thisIndex[7];
628    regionValue0 = region[iRow0];
629    regionValue1 = region[iRow1];
630    region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
631    region[iRow1] = regionValue1 - thisElement[7] * pivotValue;
632    break;
633  default:
634    if ((number & 1) != 0) {
635      number--;
636      CoinSimplexInt iRow = thisIndex[number];
637      CoinFactorizationDouble regionValue = region[iRow];
638      CoinFactorizationDouble value = thisElement[number];
639      region[iRow] = regionValue - value * pivotValue;
640    }
641    for (CoinBigIndex j = number - 1; j >= 0; j -= 2) {
642      CoinSimplexInt iRow0 = thisIndex[j];
643      CoinSimplexInt iRow1 = thisIndex[j - 1];
644      CoinFactorizationDouble regionValue0 = region[iRow0];
645      CoinFactorizationDouble regionValue1 = region[iRow1];
646      region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
647      region[iRow1] = regionValue1 - thisElement[j - 1] * pivotValue;
648    }
649    break;
650  }
651#endif
652}
653inline CoinFactorizationDouble gatherUpdate(CoinSimplexInt number,
654  const CoinFactorizationDouble *COIN_RESTRICT thisElement,
655  const CoinSimplexInt *COIN_RESTRICT thisIndex,
656  CoinFactorizationDouble *COIN_RESTRICT region)
657{
658  CoinFactorizationDouble pivotValue = 0.0;
659  for (CoinBigIndex j = 0; j < number; j++) {
660    CoinFactorizationDouble value = thisElement[j];
661    CoinSimplexInt jRow = thisIndex[j];
662    value *= region[jRow];
663    pivotValue -= value;
664  }
665  return pivotValue;
666}
667// Updates part of column (FTRANL) when densish
668void CoinAbcTypeFactorization::updateColumnLDensish(CoinIndexedVector *regionSparse) const
669{
670  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
671  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
672  CoinSimplexInt number = regionSparse->getNumElements();
673#if ABC_SMALL < 3
674  CoinSimplexInt numberNonZero = 0;
675#endif
676
677  const CoinBigIndex *COIN_RESTRICT startColumn = startColumnLAddress_;
678  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowLAddress_;
679  const CoinFactorizationDouble *COIN_RESTRICT element = elementLAddress_;
680  const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
681  const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
682  CoinSimplexInt last = numberRows_;
683  assert(last == baseL_ + numberL_);
684#if ABC_DENSE_CODE > 0
685  //can take out last bit of sparse L as empty
686  last -= numberDense_;
687#endif
688  CoinSimplexInt smallestIndex = numberRowsExtra_;
689  // do easy ones
690  for (CoinSimplexInt k = 0; k < number; k++) {
691    CoinSimplexInt iPivot = regionIndex[k];
692#ifndef ABC_ORDERED_FACTORIZATION
693    CoinSimplexInt jPivot = pivotLBackwardOrder[iPivot];
694#else
695    CoinSimplexInt jPivot = iPivot;
696#endif
697#if ABC_SMALL < 3
698    if (jPivot >= baseL_)
699      smallestIndex = CoinMin(jPivot, smallestIndex);
700    else
701      regionIndex[numberNonZero++] = iPivot;
702#else
703    if (jPivot >= baseL_)
704      smallestIndex = CoinMin(jPivot, smallestIndex);
705#endif
706  }
707  instrument_start("CoinAbcFactorizationUpdateLDensish", number + (last - smallestIndex));
708  // now others
709  for (CoinSimplexInt k = smallestIndex; k < last; k++) {
710#if 0
711    for (int j=0;j<numberRows_;j+=10) {
712      for (int jj=j;jj<CoinMin(j+10,numberRows_);jj++)
713        printf("%g ",region[jj]);
714      printf("\n");
715    }
716#endif
717#ifndef ABC_ORDERED_FACTORIZATION
718    CoinSimplexInt i = pivotLOrder[k];
719#else
720    CoinSimplexInt i = k;
721#endif
722    CoinExponent expValue = ABC_EXPONENT(region[i]);
723    if (expValue) {
724      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
725        CoinBigIndex start = startColumn[k];
726        CoinBigIndex end = startColumn[k + 1];
727        instrument_add(end - start);
728        if (TEST_INT_NONZERO(end - start)) {
729          CoinFactorizationDouble pivotValue = region[i];
730#ifndef INLINE_IT
731          for (CoinBigIndex j = start; j < end; j++) {
732            CoinSimplexInt iRow = indexRow[j];
733            CoinFactorizationDouble result = region[iRow];
734            CoinFactorizationDouble value = element[j];
735            region[iRow] = result - value * pivotValue;
736          }
737#else
738          CoinAbcScatterUpdate(end - start, pivotValue, element + start, indexRow + start, region);
739#endif
740        }
741#if ABC_SMALL < 3
742        regionIndex[numberNonZero++] = i;
743      } else {
744        region[i] = 0.0;
745#endif
746      }
747    }
748  }
749  // and dense
750#if ABC_SMALL < 3
751  for (CoinSimplexInt k = last; k < numberRows_; k++) {
752#ifndef ABC_ORDERED_FACTORIZATION
753    CoinSimplexInt i = pivotLOrder[k];
754#else
755    CoinSimplexInt i = k;
756#endif
757    CoinExponent expValue = ABC_EXPONENT(region[i]);
758    if (expValue) {
759      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
760        regionIndex[numberNonZero++] = i;
761      } else {
762        region[i] = 0.0;
763      }
764    }
765  }
766  regionSparse->setNumElements(numberNonZero);
767#endif
768  instrument_end();
769}
770// Updates part of column (FTRANL) when dense
771void CoinAbcTypeFactorization::updateColumnLDense(CoinIndexedVector *regionSparse) const
772{
773  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
774  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
775  CoinSimplexInt number = regionSparse->getNumElements();
776#if ABC_SMALL < 3
777  CoinSimplexInt numberNonZero = 0;
778#endif
779
780  const CoinBigIndex *COIN_RESTRICT startColumn = startColumnLAddress_;
781  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowLAddress_;
782  const CoinFactorizationDouble *COIN_RESTRICT element = elementLAddress_;
783  const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
784  const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
785  CoinSimplexInt last = numberRows_;
786  assert(last == baseL_ + numberL_);
787#if ABC_DENSE_CODE > 0
788  //can take out last bit of sparse L as empty
789  last -= numberDense_;
790#endif
791  CoinSimplexInt smallestIndex = numberRowsExtra_;
792  // do easy ones
793  for (CoinSimplexInt k = 0; k < number; k++) {
794    CoinSimplexInt iPivot = regionIndex[k];
795#ifndef ABC_ORDERED_FACTORIZATION
796    CoinSimplexInt jPivot = pivotLBackwardOrder[iPivot];
797#else
798    CoinSimplexInt jPivot = iPivot;
799#endif
800#if ABC_SMALL < 3
801    if (jPivot >= baseL_)
802      smallestIndex = CoinMin(jPivot, smallestIndex);
803    else
804      regionIndex[numberNonZero++] = iPivot;
805#else
806    if (jPivot >= baseL_)
807      smallestIndex = CoinMin(jPivot, smallestIndex);
808#endif
809  }
810  instrument_start("CoinAbcFactorizationUpdateLDensish", number + (last - smallestIndex));
811  // now others
812  for (CoinSimplexInt k = smallestIndex; k < last; k++) {
813#ifndef ABC_ORDERED_FACTORIZATION
814    CoinSimplexInt i = pivotLOrder[k];
815#else
816    CoinSimplexInt i = k;
817#endif
818    CoinExponent expValue = ABC_EXPONENT(region[i]);
819    if (expValue) {
820      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
821        CoinBigIndex start = startColumn[k];
822        CoinBigIndex end = startColumn[k + 1];
823        instrument_add(end - start);
824        if (TEST_INT_NONZERO(end - start)) {
825          CoinFactorizationDouble pivotValue = region[i];
826#ifndef INLINE_IT
827          for (CoinBigIndex j = start; j < end; j++) {
828            CoinSimplexInt iRow = indexRow[j];
829            CoinFactorizationDouble result = region[iRow];
830            CoinFactorizationDouble value = element[j];
831            region[iRow] = result - value * pivotValue;
832          }
833#else
834          CoinAbcScatterUpdate(end - start, pivotValue, element + start, indexRow + start, region);
835#endif
836        }
837#if ABC_SMALL < 3
838        regionIndex[numberNonZero++] = i;
839      } else {
840        region[i] = 0.0;
841#endif
842      }
843    }
844  }
845  // and dense
846#if ABC_SMALL < 3
847  for (CoinSimplexInt k = last; k < numberRows_; k++) {
848#ifndef ABC_ORDERED_FACTORIZATION
849    CoinSimplexInt i = pivotLOrder[k];
850#else
851    CoinSimplexInt i = k;
852#endif
853    CoinExponent expValue = ABC_EXPONENT(region[i]);
854    if (expValue) {
855      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
856        regionIndex[numberNonZero++] = i;
857      } else {
858        region[i] = 0.0;
859      }
860    }
861  }
862  regionSparse->setNumElements(numberNonZero);
863#endif
864  instrument_end();
865}
866#if ABC_SMALL < 2
867// Updates part of column (FTRANL) when sparse
868void CoinAbcTypeFactorization::updateColumnLSparse(CoinIndexedVector *regionSparse
869#if ABC_PARALLEL
870  ,
871  int whichSparse
872#endif
873  ) const
874{
875  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
876  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
877  CoinSimplexInt number = regionSparse->getNumElements();
878  CoinSimplexInt numberNonZero;
879
880  numberNonZero = 0;
881
882  const CoinBigIndex *COIN_RESTRICT startColumn = startColumnLAddress_;
883  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowLAddress_;
884  const CoinFactorizationDouble *COIN_RESTRICT element = elementLAddress_;
885  const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
886  CoinSimplexInt nList;
887  // use sparse_ as temporary area
888  // need to redo if this fails (just means using CoinAbcStack to compute sizes)
889  assert(sizeof(CoinSimplexInt) == sizeof(CoinBigIndex));
890  CoinAbcStack *COIN_RESTRICT stackList = reinterpret_cast< CoinAbcStack * >(sparseAddress_);
891  CoinSimplexInt *COIN_RESTRICT list = listAddress_;
892#define DO_CHAR1 1
893#if DO_CHAR1 // CHAR
894  CoinCheckZero *COIN_RESTRICT mark = markRowAddress_;
895#else
896  //BIT
897  CoinSimplexUnsignedInt *COIN_RESTRICT mark = reinterpret_cast< CoinSimplexUnsignedInt * >(markRowAddress_);
898#endif
899#if ABC_PARALLEL
900  //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
901  if (whichSparse) {
902    //printf("alternative sparse\n");
903    int addAmount = whichSparse * sizeSparseArray_;
904    stackList = reinterpret_cast< CoinAbcStack * >(reinterpret_cast< char * >(stackList) + addAmount);
905    list = reinterpret_cast< CoinSimplexInt * >(reinterpret_cast< char * >(list) + addAmount);
906#if DO_CHAR1 // CHAR
907    mark = reinterpret_cast< CoinCheckZero * >(reinterpret_cast< char * >(mark) + addAmount);
908#else
909    mark = reinterpret_cast< CoinSimplexUnsignedInt * >(reinterpret_cast< char * >(mark) + addAmount);
910#endif
911  }
912#endif
913  // mark known to be zero
914#ifdef COIN_DEBUG
915#if DO_CHAR1 // CHAR
916  for (CoinSimplexInt i = 0; i < maximumRowsExtra_; i++) {
917    assert(!mark[i]);
918  }
919#else
920  //BIT
921  for (CoinSimplexInt i = 0; i< numberRows_ >> COINFACTORIZATION_SHIFT_PER_INT; i++) {
922    assert(!mark[i]);
923  }
924#endif
925#endif
926  nList = 0;
927  for (CoinSimplexInt k = 0; k < number; k++) {
928    CoinSimplexInt kPivot = regionIndex[k];
929#ifndef ABC_ORDERED_FACTORIZATION
930    CoinSimplexInt iPivot = pivotLBackwardOrder[kPivot];
931#else
932    CoinSimplexInt iPivot = kPivot;
933#endif
934    if (iPivot >= baseL_) {
935#if DO_CHAR1 // CHAR
936      CoinCheckZero mark_k = mark[kPivot];
937#else
938      CoinSimplexUnsignedInt wordK = kPivot >> COINFACTORIZATION_SHIFT_PER_INT;
939      CoinSimplexUnsignedInt bitK = (kPivot & COINFACTORIZATION_MASK_PER_INT);
940      CoinSimplexUnsignedInt mark_k = (mark[wordK] >> bitK) & 1;
941#endif
942      if (!mark_k) {
943        CoinBigIndex j = startColumn[iPivot + 1] - 1;
944        stackList[0].stack = kPivot;
945        CoinBigIndex start = startColumn[iPivot];
946        stackList[0].next = j;
947        stackList[0].start = start;
948        CoinSimplexInt nStack = 0;
949        while (nStack >= 0) {
950          /* take off stack */
951#ifndef ABC_ORDERED_FACTORIZATION
952          iPivot = pivotLBackwardOrder[kPivot]; // put startColumn on stackstuff
953#else
954          iPivot = kPivot;
955#endif
956#if 1 //0 // what went wrong??
957          CoinBigIndex startThis = startColumn[iPivot];
958          for (; j >= startThis; j--) {
959            CoinSimplexInt jPivot = indexRow[j];
960#if DO_CHAR1 // CHAR
961            CoinCheckZero mark_j = mark[jPivot];
962#else
963            CoinSimplexUnsignedInt word0 = jPivot >> COINFACTORIZATION_SHIFT_PER_INT;
964            CoinSimplexUnsignedInt bit0 = (jPivot & COINFACTORIZATION_MASK_PER_INT);
965            CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 1;
966#endif
967            if (!mark_j) {
968#if DO_CHAR1 // CHAR
969              mark[jPivot] = 1;
970#else
971              mark[word0] |= (1 << bit0);
972#endif
973              /* and new one */
974              kPivot = jPivot;
975              break;
976            }
977          }
978          if (j >= startThis) {
979            /* put back on stack */
980            stackList[nStack].next = j - 1;
981#ifndef ABC_ORDERED_FACTORIZATION
982            iPivot = pivotLBackwardOrder[kPivot];
983#else
984            iPivot = kPivot;
985#endif
986            j = startColumn[iPivot + 1] - 1;
987            nStack++;
988            stackList[nStack].stack = kPivot;
989            assert(kPivot < numberRowsExtra_);
990            CoinBigIndex start = startColumn[iPivot];
991            stackList[nStack].next = j;
992            stackList[nStack].start = start;
993#else
994          if (j >= startColumn[iPivot] /*stackList[nStack].start*/) {
995            CoinSimplexInt jPivot = indexRow[j--];
996            /* put back on stack */
997            stackList[nStack].next = j;
998#if DO_CHAR1 // CHAR
999            CoinCheckZero mark_j = mark[jPivot];
1000#else
1001            CoinSimplexUnsignedInt word0 = jPivot >> COINFACTORIZATION_SHIFT_PER_INT;
1002            CoinSimplexUnsignedInt bit0 = (jPivot & COINFACTORIZATION_MASK_PER_INT);
1003            CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 1;
1004#endif
1005            if (!mark_j) {
1006              /* and new one */
1007              kPivot = jPivot;
1008#ifndef ABC_ORDERED_FACTORIZATION
1009              iPivot = pivotLBackwardOrder[kPivot];
1010#else
1011              iPivot = kPivot;
1012#endif
1013              j = startColumn[iPivot + 1] - 1;
1014              nStack++;
1015              stackList[nStack].stack = kPivot;
1016              assert(kPivot < numberRowsExtra_);
1017              CoinBigIndex start = startColumn[iPivot];
1018              stackList[nStack].next = j;
1019              stackList[nStack].start = start;
1020#if DO_CHAR1 // CHAR
1021              mark[jPivot] = 1;
1022#else
1023              mark[word0] |= (1 << bit0);
1024#endif
1025            }
1026#endif
1027          } else {
1028            /* finished so mark */
1029            list[nList++] = kPivot;
1030#if DO_CHAR1 // CHAR
1031            mark[kPivot] = 1;
1032#else
1033            CoinSimplexUnsignedInt wordB = kPivot >> COINFACTORIZATION_SHIFT_PER_INT;
1034            CoinSimplexUnsignedInt bitB = (kPivot & COINFACTORIZATION_MASK_PER_INT);
1035            mark[wordB] |= (1 << bitB);
1036#endif
1037            --nStack;
1038            if (nStack >= 0) {
1039              kPivot = stackList[nStack].stack;
1040#ifndef ABC_ORDERED_FACTORIZATION
1041              iPivot = pivotLBackwardOrder[kPivot];
1042#else
1043              iPivot = kPivot;
1044#endif
1045              assert(kPivot < numberRowsExtra_);
1046              j = stackList[nStack].next;
1047            }
1048          }
1049        }
1050      }
1051    } else {
1052      if (!TEST_LESS_THAN_TOLERANCE(region[kPivot])) {
1053        // just put on list
1054        regionIndex[numberNonZero++] = kPivot;
1055      } else {
1056        region[kPivot] = 0.0;
1057      }
1058    }
1059  }
1060#if 0
1061  CoinSimplexDouble temp[20000];
1062  {
1063    memcpy(temp,region,numberRows_*sizeof(CoinSimplexDouble));
1064    for (CoinSimplexInt k = 0; k < numberRows_; k++ ) {
1065      CoinSimplexInt i=pivotLOrder[k];
1066      CoinFactorizationDouble pivotValue = temp[i];
1067   
1068      if ( !TEST_LESS_THAN_TOLERANCE(pivotValue) ) {
1069        CoinBigIndex start = startColumn[k];
1070        CoinBigIndex end = startColumn[k + 1];
1071        for (CoinBigIndex j = start; j < end; j ++ ) {
1072          CoinSimplexInt iRow = indexRow[j];
1073          CoinFactorizationDouble result = temp[iRow];
1074          CoinFactorizationDouble value = element[j];
1075         
1076          temp[iRow] = result - value * pivotValue;
1077        }     
1078      } else {
1079        temp[i] = 0.0;
1080      }       
1081    }     
1082  }
1083#endif
1084  instrument_start("CoinAbcFactorizationUpdateLSparse", number + nList);
1085#if 0 //ndef NDEBUG
1086  char * test = new char[numberRows_];
1087  memset(test,0,numberRows_);
1088  for (int i=0;i<numberRows_;i++) {
1089    if (region[i])
1090      test[i]=-1;
1091  }
1092  for (CoinSimplexInt i=nList-1;i>=0;i--) {
1093    CoinSimplexInt iPivot = list[i];
1094    test[iPivot]=1;
1095    CoinSimplexInt kPivot=pivotLBackwardOrder[iPivot];
1096    CoinBigIndex start=startColumn[kPivot];
1097    CoinBigIndex end=startColumn[kPivot+1];
1098    for (CoinBigIndex j = start;j < end; j ++ ) {
1099      CoinSimplexInt iRow = indexRow[j];
1100      assert (test[iRow]<=0);
1101    }
1102  }
1103  delete [] test;
1104#endif
1105  for (CoinSimplexInt i = nList - 1; i >= 0; i--) {
1106    CoinSimplexInt iPivot = list[i];
1107#ifndef ABC_ORDERED_FACTORIZATION
1108    CoinSimplexInt kPivot = pivotLBackwardOrder[iPivot];
1109#else
1110    CoinSimplexInt kPivot = iPivot;
1111#endif
1112#if DO_CHAR1 // CHAR
1113    mark[iPivot] = 0;
1114#else
1115    CoinSimplexUnsignedInt word0 = iPivot >> COINFACTORIZATION_SHIFT_PER_INT;
1116    //CoinSimplexUnsignedInt bit0 = (iPivot & COINFACTORIZATION_MASK_PER_INT);
1117    //mark[word0]&=(~(1<<bit0));
1118    mark[word0] = 0;
1119#endif
1120    if (!TEST_LESS_THAN_TOLERANCE(region[iPivot])) {
1121      CoinFactorizationDouble pivotValue = region[iPivot];
1122      regionIndex[numberNonZero++] = iPivot;
1123      CoinBigIndex start = startColumn[kPivot];
1124      CoinBigIndex end = startColumn[kPivot + 1];
1125      instrument_add(end - start);
1126      if (TEST_INT_NONZERO(end - start)) {
1127#ifndef INLINE_IT
1128        for (CoinBigIndex j = start; j < end; j++) {
1129          CoinSimplexInt iRow = indexRow[j];
1130          CoinFactorizationDouble value = element[j];
1131          region[iRow] -= value * pivotValue;
1132        }
1133#else
1134        CoinAbcScatterUpdate(end - start, pivotValue, element + start, indexRow + start, region);
1135#endif
1136      }
1137    } else {
1138      region[iPivot] = 0.0;
1139    }
1140  }
1141#if 0
1142  {
1143    for (CoinSimplexInt k = 0; k < numberRows_; k++ ) {
1144      assert (fabs(region[k]-temp[k])<1.0e-1);
1145    }
1146  }
1147#endif
1148  regionSparse->setNumElements(numberNonZero);
1149  instrument_end_and_adjust(1.3);
1150}
1151#endif
1152#if FACTORIZATION_STATISTICS
1153extern double twoThresholdX;
1154#endif
1155CoinSimplexInt
1156CoinAbcTypeFactorization::updateTwoColumnsFT(CoinIndexedVector &regionFTX,
1157  CoinIndexedVector &regionOtherX)
1158{
1159  CoinIndexedVector *regionFT = &regionFTX;
1160  CoinIndexedVector *regionOther = &regionOtherX;
1161#if ABC_DEBUG
1162  regionFT->checkClean();
1163  regionOther->checkClean();
1164#endif
1165  toLongArray(regionFT, 0);
1166  toLongArray(regionOther, 1);
1167#ifdef ABC_ORDERED_FACTORIZATION
1168  // Permute in for Ftran
1169  permuteInForFtran(regionFTX);
1170  permuteInForFtran(regionOtherX);
1171#endif
1172  CoinSimplexInt numberNonZero = regionOther->getNumElements();
1173  CoinSimplexInt numberNonZeroFT = regionFT->getNumElements();
1174  //  ******* L
1175  updateColumnL(regionFT
1176#if ABC_SMALL < 2
1177    ,
1178    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
1179#endif
1180  );
1181  updateColumnL(regionOther
1182#if ABC_SMALL < 2
1183    ,
1184    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
1185#endif
1186  );
1187  //row bits here
1188  updateColumnR(regionFT
1189#if ABC_SMALL < 2
1190    ,
1191    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
1192#endif
1193  );
1194  updateColumnR(regionOther
1195#if ABC_SMALL < 2
1196    ,
1197    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
1198#endif
1199  );
1200  bool doFT = storeFT(regionFT);
1201  //update counts
1202  //  ******* U - see if densish
1203#if ABC_SMALL < 2
1204  // Guess at number at end
1205  CoinSimplexInt goSparse = 0;
1206  if (gotSparse()) {
1207    CoinSimplexInt numberNonZero = (regionOther->getNumElements() + regionFT->getNumElements()) >> 1;
1208    double average = 0.25 * (ftranAverageAfterL_ * twiddleFtranFactor1() + ftranFTAverageAfterL_ * twiddleFtranFTFactor1());
1209    assert(average);
1210    CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * average);
1211    if (newNumber < sparseThreshold_)
1212      goSparse = 2;
1213  }
1214#if FACTORIZATION_STATISTICS
1215  int twoThreshold = twoThresholdX;
1216#else
1217#define twoThreshold 1000
1218#endif
1219#else
1220#define goSparse false
1221#define twoThreshold 1000
1222#endif
1223  if (!goSparse && numberRows_ < twoThreshold) {
1224    CoinFactorizationDouble *COIN_RESTRICT arrayFT = denseVector(regionFT);
1225    CoinSimplexInt *COIN_RESTRICT indexFT = regionFT->getIndices();
1226    CoinFactorizationDouble *COIN_RESTRICT arrayUpdate = denseVector(regionOther);
1227    CoinSimplexInt *COIN_RESTRICT indexUpdate = regionOther->getIndices();
1228    updateTwoColumnsUDensish(numberNonZeroFT, arrayFT, indexFT,
1229      numberNonZero, arrayUpdate, indexUpdate);
1230    regionFT->setNumElements(numberNonZeroFT);
1231    regionOther->setNumElements(numberNonZero);
1232  } else {
1233    // sparse
1234    updateColumnU(regionFT
1235#if ABC_SMALL < 2
1236      ,
1237      reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
1238#endif
1239    );
1240    numberNonZeroFT = regionFT->getNumElements();
1241    updateColumnU(regionOther
1242#if ABC_SMALL < 2
1243      ,
1244      reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
1245#endif
1246    );
1247    numberNonZero = regionOther->getNumElements();
1248  }
1249  fromLongArray(0);
1250  fromLongArray(1);
1251#if ABC_DEBUG
1252  regionFT->checkClean();
1253  regionOther->checkClean();
1254#endif
1255  if (doFT)
1256    return numberNonZeroFT;
1257  else
1258    return -numberNonZeroFT;
1259}
1260// Updates part of 2 columns (FTRANU) real work
1261void CoinAbcTypeFactorization::updateTwoColumnsUDensish(
1262  CoinSimplexInt &numberNonZero1,
1263  CoinFactorizationDouble *COIN_RESTRICT region1,
1264  CoinSimplexInt *COIN_RESTRICT index1,
1265  CoinSimplexInt &numberNonZero2,
1266  CoinFactorizationDouble *COIN_RESTRICT region2,
1267  CoinSimplexInt *COIN_RESTRICT index2) const
1268{
1269#ifdef ABC_USE_FUNCTION_POINTERS
1270  scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1271  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1272#else
1273  const CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1274  const CoinBigIndex *COIN_RESTRICT startColumn = startColumnUAddress_;
1275  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1276  const CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1277#endif
1278  CoinSimplexInt numberNonZeroA = 0;
1279  CoinSimplexInt numberNonZeroB = 0;
1280  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1281  const CoinSimplexInt *COIN_RESTRICT pivotLinked = pivotLinkedBackwardsAddress_;
1282  CoinSimplexInt jRow = pivotLinked[numberRows_];
1283  instrument_start("CoinAbcFactorizationUpdateTwoUDense", 2 * numberRows_);
1284#if 1 //def DONT_USE_SLACKS
1285  while (jRow != -1 /*lastSlack_*/) {
1286#else
1287  // would need extra code
1288  while (jRow != lastSlack_) {
1289#endif
1290    bool nonZero2 = (TEST_DOUBLE_NONZERO(region2[jRow]));
1291    bool nonZero1 = (TEST_DOUBLE_NONZERO(region1[jRow]));
1292#ifndef ABC_USE_FUNCTION_POINTERS
1293    int numberIn = numberInColumn[jRow];
1294#else
1295    scatterStruct &COIN_RESTRICT scatter = scatterPointer[jRow];
1296    CoinSimplexInt numberIn = scatter.number;
1297#endif
1298    if (nonZero1 || nonZero2) {
1299#ifndef ABC_USE_FUNCTION_POINTERS
1300      CoinBigIndex start = startColumn[jRow];
1301      const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
1302      const CoinSimplexInt *COIN_RESTRICT thisIndex = indexRow + start;
1303#else
1304      const CoinFactorizationDouble *COIN_RESTRICT thisElement = elementUColumnPlus + scatter.offset;
1305      const CoinSimplexInt *COIN_RESTRICT thisIndex = reinterpret_cast< const CoinSimplexInt * >(thisElement + numberIn);
1306#endif
1307      CoinFactorizationDouble pivotValue2 = region2[jRow];
1308      CoinFactorizationDouble pivotMult = pivotRegion[jRow];
1309      assert(pivotMult);
1310      CoinFactorizationDouble pivotValue2a = pivotValue2 * pivotMult;
1311      bool do2 = !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2);
1312      region2[jRow] = 0.0;
1313      CoinFactorizationDouble pivotValue1 = region1[jRow];
1314      region1[jRow] = 0.0;
1315      CoinFactorizationDouble pivotValue1a = pivotValue1 * pivotMult;
1316      bool do1 = !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue1);
1317      if (do2) {
1318        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2a)) {
1319          region2[jRow] = pivotValue2a;
1320          index2[numberNonZeroB++] = jRow;
1321        }
1322      }
1323      if (do1) {
1324        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue1a)) {
1325          region1[jRow] = pivotValue1a;
1326          index1[numberNonZeroA++] = jRow;
1327        }
1328      }
1329      instrument_add(numberIn);
1330      if (numberIn) {
1331        if (do2) {
1332          if (!do1) {
1333            // just region 2
1334            for (CoinBigIndex j = numberIn - 1; j >= 0; j--) {
1335              CoinSimplexInt iRow = thisIndex[j];
1336              CoinFactorizationDouble value = thisElement[j];
1337              assert(value);
1338#ifdef NO_LOAD
1339              region2[iRow] -= value * pivotValue2;
1340#else
1341              CoinFactorizationDouble regionValue2 = region2[iRow];
1342              region2[iRow] = regionValue2 - value * pivotValue2;
1343#endif
1344            }
1345          } else {
1346            // both
1347            instrument_add(numberIn);
1348            for (CoinBigIndex j = numberIn - 1; j >= 0; j--) {
1349              CoinSimplexInt iRow = thisIndex[j];
1350              CoinFactorizationDouble value = thisElement[j];
1351#ifdef NO_LOAD
1352              region1[iRow] -= value * pivotValue1;
1353              region2[iRow] -= value * pivotValue2;
1354#else
1355              CoinFactorizationDouble regionValue1 = region1[iRow];
1356              CoinFactorizationDouble regionValue2 = region2[iRow];
1357              region1[iRow] = regionValue1 - value * pivotValue1;
1358              region2[iRow] = regionValue2 - value * pivotValue2;
1359#endif
1360            }
1361          }
1362        } else if (do1) {
1363          // just region 1
1364          for (CoinBigIndex j = numberIn - 1; j >= 0; j--) {
1365            CoinSimplexInt iRow = thisIndex[j];
1366            CoinFactorizationDouble value = thisElement[j];
1367            assert(value);
1368#ifdef NO_LOAD
1369            region1[iRow] -= value * pivotValue1;
1370#else
1371            CoinFactorizationDouble regionValue1 = region1[iRow];
1372            region1[iRow] = regionValue1 - value * pivotValue1;
1373#endif
1374          }
1375        }
1376      } else {
1377        // no elements in column
1378      }
1379    }
1380    jRow = pivotLinked[jRow];
1381  }
1382  numberNonZero1 = numberNonZeroA;
1383  numberNonZero2 = numberNonZeroB;
1384#if ABC_SMALL < 2
1385  if (factorizationStatistics()) {
1386    ftranFTCountAfterU_ += numberNonZeroA;
1387    ftranCountAfterU_ += numberNonZeroB;
1388  }
1389#endif
1390  instrument_end();
1391}
1392//  updateColumnU.  Updates part of column (FTRANU)
1393void CoinAbcTypeFactorization::updateColumnU(CoinIndexedVector *regionSparse
1394#if ABC_SMALL < 2
1395  ,
1396  CoinAbcStatistics &statistics
1397#endif
1398#if ABC_PARALLEL
1399  ,
1400  int whichSparse
1401#endif
1402  ) const
1403{
1404#if CILK_CONFLICT > 0
1405#if ABC_PARALLEL
1406  // for conflicts
1407#if CILK_CONFLICT > 1
1408  printf("file %s line %d which %d\n", __FILE__, __LINE__, whichSparse);
1409#endif
1410  abc_assert((cilk_conflict & (1 << whichSparse)) == 0);
1411  cilk_conflict |= (1 << whichSparse);
1412#else
1413  abc_assert((cilk_conflict & (1 << 0)) == 0);
1414  cilk_conflict |= (1 << 0);
1415#endif
1416#endif
1417#if ABC_SMALL < 2
1418  CoinSimplexInt numberNonZero = regionSparse->getNumElements();
1419  int goSparse;
1420  // Guess at number at end
1421  if (gotSparse()) {
1422    double average = statistics.averageAfterU_ * twiddleFactor1S();
1423    assert(average);
1424    CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * average);
1425    if (newNumber < sparseThreshold_)
1426      goSparse = 1;
1427    else if (numberNonZero * 3 < numberRows_)
1428      goSparse = 0;
1429    else
1430      goSparse = -1;
1431  } else {
1432    goSparse = 0;
1433  }
1434#endif
1435  if (!goSparse) {
1436    // densish
1437    updateColumnUDensish(regionSparse);
1438#if ABC_SMALL < 2
1439  } else if (goSparse < 0) {
1440    // dense
1441    updateColumnUDense(regionSparse);
1442  } else {
1443    // sparse
1444    updateColumnUSparse(regionSparse
1445#if ABC_PARALLEL
1446      ,
1447      whichSparse
1448#endif
1449    );
1450#endif
1451  }
1452#if ABC_SMALL < 2
1453  if (factorizationStatistics()) {
1454    statistics.countAfterU_ += regionSparse->getNumElements();
1455  }
1456#endif
1457#if CILK_CONFLICT > 0
1458#if ABC_PARALLEL
1459  // for conflicts
1460  abc_assert((cilk_conflict & (1 << whichSparse)) != 0);
1461  cilk_conflict &= ~(1 << whichSparse);
1462#else
1463  abc_assert((cilk_conflict & (1 << 0)) != 0);
1464  cilk_conflict &= ~(1 << 0);
1465#endif
1466#endif
1467}
1468//#define COUNT_U
1469#ifdef COUNT_U
1470static int nUDense[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
1471static int nUSparse[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
1472#endif
1473//  updateColumnU.  Updates part of column (FTRANU)
1474// Updates part of column (FTRANU) real work
1475void CoinAbcTypeFactorization::updateColumnUDensish(CoinIndexedVector *regionSparse) const
1476{
1477  instrument_start("CoinAbcFactorizationUpdateUDensish", 2 * numberRows_);
1478  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1479  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
1480  CoinSimplexInt numberNonZero = 0;
1481  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1482#ifdef ABC_USE_FUNCTION_POINTERS
1483  scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1484  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1485#else
1486  const CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1487  const CoinBigIndex *COIN_RESTRICT startColumn = startColumnUAddress_;
1488  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1489  const CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1490#endif
1491
1492  const CoinSimplexInt *COIN_RESTRICT pivotLinked = pivotLinkedBackwardsAddress_;
1493  CoinSimplexInt jRow = pivotLinked[numberRows_];
1494#define ETAS_1 2
1495#define TEST_BEFORE
1496#ifdef TEST_BEFORE
1497  CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1498#endif
1499#ifdef DONT_USE_SLACKS
1500  while (jRow != -1 /*lastSlack_*/) {
1501#else
1502  while (jRow != lastSlack_) {
1503#endif
1504#ifndef TEST_BEFORE
1505    CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1506#endif
1507    if (expValue) {
1508      if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
1509        CoinFactorizationDouble pivotValue = region[jRow];
1510#if ETAS_1 > 1
1511        CoinFactorizationDouble pivotValue2 = pivotValue * pivotRegion[jRow];
1512#endif
1513#ifndef ABC_USE_FUNCTION_POINTERS
1514        int number = numberInColumn[jRow];
1515#else
1516        scatterStruct &COIN_RESTRICT scatter = scatterPointer[jRow];
1517        CoinSimplexInt number = scatter.number;
1518#endif
1519        instrument_add(number);
1520        if (TEST_INT_NONZERO(number)) {
1521#ifdef COUNT_U
1522          {
1523            int k = numberInColumn[jRow];
1524            if (k > 10)
1525              k = 11;
1526            nUDense[k]++;
1527            int kk = 0;
1528            for (int k = 0; k < 12; k++)
1529              kk += nUDense[k];
1530            if (kk % 1000000 == 0) {
1531              printf("ZZ");
1532              for (int k = 0; k < 12; k++)
1533                printf(" %d", nUDense[k]);
1534              printf("\n");
1535            }
1536          }
1537#endif
1538#ifndef ABC_USE_FUNCTION_POINTERS
1539          CoinBigIndex start = startColumn[jRow];
1540#ifndef INLINE_IT
1541          const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
1542          const CoinSimplexInt *COIN_RESTRICT thisIndex = indexRow + start;
1543          for (CoinBigIndex j = number - 1; j >= 0; j--) {
1544            CoinSimplexInt iRow = thisIndex[j];
1545            CoinFactorizationDouble regionValue = region[iRow];
1546            CoinFactorizationDouble value = thisElement[j];
1547            assert(value);
1548            region[iRow] = regionValue - value * pivotValue;
1549          }
1550#else
1551          CoinAbcScatterUpdate(number, pivotValue, element + start, indexRow + start, region);
1552#endif
1553#else
1554          CoinBigIndex start = scatter.offset;
1555#if ABC_USE_FUNCTION_POINTERS
1556          (*(scatter.functionPointer))(number, pivotValue, elementUColumnPlus + start, region);
1557#else
1558          CoinAbcScatterUpdate(number, pivotValue, elementUColumnPlus + start, region);
1559#endif
1560#endif
1561        }
1562        CoinSimplexInt kRow = jRow;
1563        jRow = pivotLinked[jRow];
1564#ifdef TEST_BEFORE
1565        expValue = ABC_EXPONENT(region[jRow]);
1566#endif
1567#if ETAS_1 < 2
1568        pivotValue *= pivotRegion[kRow];
1569        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
1570          region[kRow] = pivotValue;
1571          regionIndex[numberNonZero++] = kRow;
1572        } else {
1573          region[kRow] = 0.0;
1574        }
1575#else
1576        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2)) {
1577          region[kRow] = pivotValue2;
1578          regionIndex[numberNonZero++] = kRow;
1579        } else {
1580          region[kRow] = 0.0;
1581        }
1582#endif
1583      } else {
1584        CoinSimplexInt kRow = jRow;
1585        jRow = pivotLinked[jRow];
1586#ifdef TEST_BEFORE
1587        expValue = ABC_EXPONENT(region[jRow]);
1588#endif
1589        region[kRow] = 0.0;
1590      }
1591    } else {
1592      jRow = pivotLinked[jRow];
1593#ifdef TEST_BEFORE
1594      expValue = ABC_EXPONENT(region[jRow]);
1595#endif
1596    }
1597  }
1598#ifndef DONT_USE_SLACKS
1599  while (jRow != -1) {
1600#ifndef TEST_BEFORE
1601    CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1602#endif
1603    if (expValue) {
1604      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
1605#if SLACK_VALUE == -1
1606        CoinFactorizationDouble pivotValue = region[jRow];
1607        assert(pivotRegion[jRow] == SLACK_VALUE);
1608        region[jRow] = -pivotValue;
1609#endif
1610        regionIndex[numberNonZero++] = jRow;
1611      } else {
1612        region[jRow] = 0.0;
1613      }
1614    }
1615    jRow = pivotLinked[jRow];
1616#ifdef TEST_BEFORE
1617    expValue = ABC_EXPONENT(region[jRow]);
1618#endif
1619  }
1620#endif
1621  regionSparse->setNumElements(numberNonZero);
1622  instrument_end();
1623}
1624//  updateColumnU.  Updates part of column (FTRANU)
1625// Updates part of column (FTRANU) real work
1626void CoinAbcTypeFactorization::updateColumnUDense(CoinIndexedVector *regionSparse) const
1627{
1628  instrument_start("CoinAbcFactorizationUpdateUDensish", 2 * numberRows_);
1629  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1630  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
1631  const CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1632  CoinSimplexInt numberNonZero = 0;
1633  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1634
1635#ifdef ABC_USE_FUNCTION_POINTERS
1636  scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1637  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1638#else
1639  const CoinBigIndex *COIN_RESTRICT startColumn = startColumnUAddress_;
1640  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1641  const CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1642#endif
1643  const CoinSimplexInt *COIN_RESTRICT pivotLinked = pivotLinkedBackwardsAddress_;
1644  CoinSimplexInt jRow = pivotLinked[numberRows_];
1645#define ETAS_1 2
1646#define TEST_BEFORE
1647#ifdef TEST_BEFORE
1648  CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1649#endif
1650  //int ixxxxxx=0;
1651#ifdef DONT_USE_SLACKS
1652  while (jRow != -1 /*lastSlack_*/) {
1653#else
1654  while (jRow != lastSlack_) {
1655#endif
1656#if 0
1657    double largest=1.0;
1658    int iLargest=-1;
1659    ixxxxxx++;
1660    for (int i=0;i<numberRows_;i++) {
1661      if (fabs(region[i])>largest) {
1662        largest=fabs(region[i]);
1663        iLargest=i;
1664      }
1665    }
1666    if (iLargest>=0)
1667      printf("largest %g on row %d after %d\n",largest,iLargest,ixxxxxx);
1668#endif
1669#ifndef TEST_BEFORE
1670    CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1671#endif
1672    if (expValue) {
1673      if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
1674        CoinFactorizationDouble pivotValue = region[jRow];
1675#if ETAS_1 > 1
1676        CoinFactorizationDouble pivotValue2 = pivotValue * pivotRegion[jRow];
1677#endif
1678#ifndef ABC_USE_FUNCTION_POINTERS
1679        int number = numberInColumn[jRow];
1680#else
1681        scatterStruct &COIN_RESTRICT scatter = scatterPointer[jRow];
1682        CoinSimplexInt number = scatter.number;
1683#endif
1684        instrument_add(number);
1685        if (TEST_INT_NONZERO(number)) {
1686#ifdef COUNT_U
1687          {
1688            int k = numberInColumn[jRow];
1689            if (k > 10)
1690              k = 11;
1691            nUDense[k]++;
1692            int kk = 0;
1693            for (int k = 0; k < 12; k++)
1694              kk += nUDense[k];
1695            if (kk % 1000000 == 0) {
1696              printf("ZZ");
1697              for (int k = 0; k < 12; k++)
1698                printf(" %d", nUDense[k]);
1699              printf("\n");
1700            }
1701          }
1702#endif
1703#ifndef ABC_USE_FUNCTION_POINTERS
1704          CoinBigIndex start = startColumn[jRow];
1705#ifndef INLINE_IT
1706          const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
1707          const CoinSimplexInt *COIN_RESTRICT thisIndex = indexRow + start;
1708          for (CoinBigIndex j = number - 1; j >= 0; j--) {
1709            CoinSimplexInt iRow = thisIndex[j];
1710            CoinFactorizationDouble regionValue = region[iRow];
1711            CoinFactorizationDouble value = thisElement[j];
1712            assert(value);
1713            region[iRow] = regionValue - value * pivotValue;
1714          }
1715#else
1716          CoinAbcScatterUpdate(number, pivotValue, element + start, indexRow + start, region);
1717#endif
1718#else
1719          CoinBigIndex start = scatter.offset;
1720#if ABC_USE_FUNCTION_POINTERS
1721          (*(scatter.functionPointer))(number, pivotValue, elementUColumnPlus + start, region);
1722#else
1723          CoinAbcScatterUpdate(number, pivotValue, elementUColumnPlus + start, region);
1724#endif
1725#endif
1726        }
1727        CoinSimplexInt kRow = jRow;
1728        jRow = pivotLinked[jRow];
1729#ifdef TEST_BEFORE
1730        expValue = ABC_EXPONENT(region[jRow]);
1731#endif
1732#if ETAS_1 < 2
1733        pivotValue *= pivotRegion[kRow];
1734        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
1735          region[kRow] = pivotValue;
1736          regionIndex[numberNonZero++] = kRow;
1737        } else {
1738          region[kRow] = 0.0;
1739        }
1740#else
1741        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2)) {
1742          region[kRow] = pivotValue2;
1743          regionIndex[numberNonZero++] = kRow;
1744        } else {
1745          region[kRow] = 0.0;
1746        }
1747#endif
1748      } else {
1749        CoinSimplexInt kRow = jRow;
1750        jRow = pivotLinked[jRow];
1751#ifdef TEST_BEFORE
1752        expValue = ABC_EXPONENT(region[jRow]);
1753#endif
1754        region[kRow] = 0.0;
1755      }
1756    } else {
1757      jRow = pivotLinked[jRow];
1758#ifdef TEST_BEFORE
1759      expValue = ABC_EXPONENT(region[jRow]);
1760#endif
1761    }
1762  }
1763#ifndef DONT_USE_SLACKS
1764  while (jRow != -1) {
1765#ifndef TEST_BEFORE
1766    CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1767#endif
1768    if (expValue) {
1769      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
1770#if SLACK_VALUE == -1
1771        CoinFactorizationDouble pivotValue = region[jRow];
1772        assert(pivotRegion[jRow] == SLACK_VALUE);
1773        region[jRow] = -pivotValue;
1774#endif
1775        regionIndex[numberNonZero++] = jRow;
1776      } else {
1777        region[jRow] = 0.0;
1778      }
1779    }
1780    jRow = pivotLinked[jRow];
1781#ifdef TEST_BEFORE
1782    expValue = ABC_EXPONENT(region[jRow]);
1783#endif
1784  }
1785#endif
1786  regionSparse->setNumElements(numberNonZero);
1787  instrument_end();
1788}
1789#if ABC_SMALL < 2
1790//  updateColumnU.  Updates part of column (FTRANU)
1791/*
1792  Since everything is in order I should be able to do a better job of
1793  marking stuff - think.  Also as L is static maybe I can do something
1794  better there (I know I could if I marked the depth of every element
1795  but that would lead to other inefficiencies.
1796*/
1797void CoinAbcTypeFactorization::updateColumnUSparse(CoinIndexedVector *regionSparse
1798#if ABC_PARALLEL
1799  ,
1800  int whichSparse
1801#endif
1802  ) const
1803{
1804  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1805  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
1806  CoinSimplexInt numberNonZero = regionSparse->getNumElements();
1807  //const CoinBigIndex * COIN_RESTRICT startColumn = startColumnUAddress_;
1808  //const CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1809  //const CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
1810  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1811  // use sparse_ as temporary area
1812  // mark known to be zero
1813#define COINFACTORIZATION_SHIFT_PER_INT2 (COINFACTORIZATION_SHIFT_PER_INT - 1)
1814#define COINFACTORIZATION_MASK_PER_INT2 (COINFACTORIZATION_MASK_PER_INT >> 1)
1815  // need to redo if this fails (just means using CoinAbcStack to compute sizes)
1816  assert(sizeof(CoinSimplexInt) == sizeof(CoinBigIndex));
1817  CoinAbcStack *COIN_RESTRICT stackList = reinterpret_cast< CoinAbcStack * >(sparseAddress_);
1818  CoinSimplexInt *COIN_RESTRICT list = listAddress_;
1819#ifndef ABC_USE_FUNCTION_POINTERS
1820  const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1821#else
1822  scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1823  const CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1824  const CoinSimplexInt *COIN_RESTRICT indexRow = reinterpret_cast< const CoinSimplexInt * >(elementUColumnPlusAddress_);
1825#endif
1826  //#define DO_CHAR2 1
1827#if DO_CHAR2 // CHAR
1828  CoinCheckZero *COIN_RESTRICT mark = markRowAddress_;
1829#else
1830  //BIT
1831  CoinSimplexUnsignedInt *COIN_RESTRICT mark = reinterpret_cast< CoinSimplexUnsignedInt * >(markRowAddress_);
1832#endif
1833#if ABC_PARALLEL
1834  //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
1835  if (whichSparse) {
1836    int addAmount = whichSparse * sizeSparseArray_;
1837    stackList = reinterpret_cast< CoinAbcStack * >(reinterpret_cast< char * >(stackList) + addAmount);
1838    list = reinterpret_cast< CoinSimplexInt * >(reinterpret_cast< char * >(list) + addAmount);
1839#if DO_CHAR2 // CHAR
1840    mark = reinterpret_cast< CoinCheckZero * >(reinterpret_cast< char * >(mark) + addAmount);
1841#else
1842    mark = reinterpret_cast< CoinSimplexUnsignedInt * >(reinterpret_cast< char * >(mark) + addAmount);
1843#endif
1844  }
1845#endif
1846#ifdef COIN_DEBUG
1847#if DO_CHAR2 // CHAR
1848  for (CoinSimplexInt i = 0; i < maximumRowsExtra_; i++) {
1849    assert(!mark[i]);
1850  }
1851#else
1852  //BIT
1853  for (CoinSimplexInt i = 0; i< numberRows_ >> COINFACTORIZATION_SHIFT_PER_INT; i++) {
1854    assert(!mark[i]);
1855  }
1856#endif
1857#endif
1858  CoinSimplexInt nList = 0;
1859  // move slacks to end of stack list
1860  int *COIN_RESTRICT putLast = list + numberRows_;
1861  int *COIN_RESTRICT put = putLast;
1862  for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
1863    CoinSimplexInt kPivot = regionIndex[i];
1864#if DO_CHAR2 // CHAR
1865    CoinCheckZero mark_B = mark[kPivot];
1866#else
1867    CoinSimplexUnsignedInt wordB = kPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1868    CoinSimplexUnsignedInt bitB = (kPivot & COINFACTORIZATION_MASK_PER_INT2) << 1;
1869    CoinSimplexUnsignedInt mark_B = (mark[wordB] >> bitB) & 3;
1870#endif
1871    if (!mark_B) {
1872#if DO_CHAR2 // CHAR
1873      mark[kPivot] = 1;
1874#else
1875      mark[wordB] |= (1 << bitB);
1876#endif
1877#ifndef ABC_USE_FUNCTION_POINTERS
1878      CoinBigIndex start = startColumn[kPivot];
1879      CoinSimplexInt number = numberInColumn[kPivot];
1880#else
1881      scatterStruct &COIN_RESTRICT scatter = scatterPointer[kPivot];
1882      CoinSimplexInt number = scatter.number;
1883      CoinBigIndex start = (scatter.offset + number) << 1;
1884#endif
1885      stackList[0].stack = kPivot;
1886      stackList[0].next = start + number;
1887      stackList[0].start = start;
1888      CoinSimplexInt nStack = 0;
1889      while (nStack >= 0) {
1890        /* take off stack */
1891        CoinBigIndex j = stackList[nStack].next - 1;
1892        CoinBigIndex start = stackList[nStack].start;
1893#if DO_CHAR2 == 0 // CHAR
1894        CoinSimplexUnsignedInt word0;
1895        CoinSimplexUnsignedInt bit0;
1896#endif
1897        CoinSimplexInt jPivot;
1898        for (; j >= start; j--) {
1899          jPivot = indexRow[j];
1900#if DO_CHAR2 // CHAR
1901          CoinCheckZero mark_j = mark[jPivot];
1902#else
1903          word0 = jPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1904          bit0 = (jPivot & COINFACTORIZATION_MASK_PER_INT2) << 1;
1905          CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 3;
1906#endif
1907          if (!mark_j)
1908            break;
1909        }
1910        if (j >= start) {
1911          /* put back on stack */
1912          stackList[nStack].next = j;
1913          /* and new one */
1914#ifndef ABC_USE_FUNCTION_POINTERS
1915          CoinBigIndex start = startColumn[jPivot];
1916          CoinSimplexInt number = numberInColumn[jPivot];
1917#else
1918          scatterStruct &COIN_RESTRICT scatter = scatterPointer[jPivot];
1919          CoinSimplexInt number = scatter.number;
1920          CoinBigIndex start = (scatter.offset + number) << 1;
1921#endif
1922          if (number) {
1923            nStack++;
1924            stackList[nStack].stack = jPivot;
1925            stackList[nStack].next = start + number;
1926            stackList[nStack].start = start;
1927#if DO_CHAR2 // CHAR
1928            mark[jPivot] = 1;
1929#else
1930            mark[word0] |= (1 << bit0);
1931#endif
1932          } else {
1933            // can do at once
1934#ifndef NDEBUG
1935#if DO_CHAR2 // CHAR
1936            CoinCheckZero mark_j = mark[jPivot];
1937#else
1938            CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 3;
1939#endif
1940            assert(mark_j != 3);
1941#endif
1942#if ABC_SMALL < 2
1943            if (!start) {
1944              // slack - put at end
1945              --put;
1946              *put = jPivot;
1947              assert(pivotRegion[jPivot] == 1.0);
1948            } else {
1949              list[nList++] = jPivot;
1950            }
1951#else
1952            list[nList++] = jPivot;
1953#endif
1954#if DO_CHAR2 // CHAR
1955            mark[jPivot] = 3;
1956#else
1957            mark[word0] |= (3 << bit0);
1958#endif
1959          }
1960        } else {
1961          /* finished so mark */
1962          CoinSimplexInt kPivot = stackList[nStack].stack;
1963#if DO_CHAR2 // CHAR
1964          CoinCheckZero mark_B = mark[kPivot];
1965#else
1966          CoinSimplexUnsignedInt wordB = kPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1967          CoinSimplexUnsignedInt bitB = (kPivot & COINFACTORIZATION_MASK_PER_INT2) << 1;
1968          CoinSimplexUnsignedInt mark_B = (mark[wordB] >> bitB) & 3;
1969#endif
1970          assert(mark_B != 3);
1971          //if (mark_B!=3) {
1972          list[nList++] = kPivot;
1973#if DO_CHAR2 // CHAR
1974          mark[kPivot] = 3;
1975#else
1976          mark[wordB] |= (3 << bitB);
1977#endif
1978          //}
1979          --nStack;
1980        }
1981      }
1982    }
1983  }
1984  instrument_start("CoinAbcFactorizationUpdateUSparse", numberNonZero + 2 * nList);
1985  numberNonZero = 0;
1986  list[-1] = -1;
1987  //assert (nList);
1988  for (CoinSimplexInt i = nList - 1; i >= 0; i--) {
1989    CoinSimplexInt iPivot = list[i];
1990    CoinExponent expValue = ABC_EXPONENT(region[iPivot]);
1991#if DO_CHAR2 // CHAR
1992    mark[iPivot] = 0;
1993#else
1994    CoinSimplexInt word0 = iPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1995    mark[word0] = 0;
1996#endif
1997    if (expValue) {
1998      if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
1999        CoinFactorizationDouble pivotValue = region[iPivot];
2000#if ETAS_1 > 1
2001        CoinFactorizationDouble pivotValue2 = pivotValue * pivotRegion[iPivot];
2002#endif
2003#ifndef ABC_USE_FUNCTION_POINTERS
2004        CoinSimplexInt number = numberInColumn[iPivot];
2005#else
2006        scatterStruct &COIN_RESTRICT scatter = scatterPointer[iPivot];
2007        CoinSimplexInt number = scatter.number;
2008#endif
2009        if (TEST_INT_NONZERO(number)) {
2010#ifdef COUNT_U
2011          {
2012            int k = numberInColumn[iPivot];
2013            if (k > 10)
2014              k = 11;
2015            nUSparse[k]++;
2016            int kk = 0;
2017            for (int k = 0; k < 12; k++)
2018              kk += nUSparse[k];
2019            if (kk % 1000000 == 0) {
2020              printf("ZZsparse");
2021              for (int k = 0; k < 12; k++)
2022                printf(" %d", nUSparse[k]);
2023              printf("\n");
2024            }
2025          }
2026#endif
2027#ifndef ABC_USE_FUNCTION_POINTERS
2028          CoinBigIndex start = startColumn[iPivot];
2029#else
2030          //CoinBigIndex start = startColumn[iPivot];
2031          CoinBigIndex start = scatter.offset;
2032#endif
2033          instrument_add(number);
2034#ifndef INLINE_IT
2035          CoinBigIndex j;
2036          for (j = start; j < start + number; j++) {
2037            CoinFactorizationDouble value = element[j];
2038            assert(value);
2039            CoinSimplexInt iRow = indexRow[j];
2040            region[iRow] -= value * pivotValue;
2041          }
2042#else
2043#ifdef ABC_USE_FUNCTION_POINTERS
2044#if ABC_USE_FUNCTION_POINTERS
2045          (*(scatter.functionPointer))(number, pivotValue, elementUColumnPlus + start, region);
2046#else
2047          CoinAbcScatterUpdate(number, pivotValue, elementUColumnPlus + start, region);
2048#endif
2049#else
2050          CoinAbcScatterUpdate(number, pivotValue, element + start, indexRow + start, region);
2051#endif
2052#endif
2053        }
2054#if ETAS_1 < 2
2055        pivotValue *= pivotRegion[iPivot];
2056        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
2057          region[iPivot] = pivotValue;
2058          regionIndex[numberNonZero++] = iPivot;
2059        }
2060#else
2061        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2)) {
2062          region[iPivot] = pivotValue2;
2063          regionIndex[numberNonZero++] = iPivot;
2064        } else {
2065          region[iPivot] = 0.0;
2066        }
2067#endif
2068      } else {
2069        region[iPivot] = 0.0;
2070      }
2071    }
2072  }
2073#if ABC_SMALL < 2
2074  // slacks
2075  for (; put < putLast; put++) {
2076    int iPivot = *put;
2077    CoinExponent expValue = ABC_EXPONENT(region[iPivot]);
2078#if DO_CHAR2 // CHAR
2079    mark[iPivot] = 0;
2080#else
2081    CoinSimplexInt word0 = iPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
2082    mark[word0] = 0;
2083#endif
2084    if (expValue) {
2085      if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
2086        CoinFactorizationDouble pivotValue = region[iPivot];
2087        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
2088          region[iPivot] = pivotValue;
2089          regionIndex[numberNonZero++] = iPivot;
2090        }
2091      } else {
2092        region[iPivot] = 0.0;
2093      }
2094    }
2095  }
2096#endif
2097#ifdef COIN_DEBUG
2098  for (CoinSimplexInt i = 0; i< maximumRowsExtra_ >> COINFACTORIZATION_SHIFT_PER_INT2; i++) {
2099    assert(!mark[i]);
2100  }
2101#endif
2102  regionSparse->setNumElements(numberNonZero);
2103  instrument_end_and_adjust(1.3);
2104}
2105#endif
2106// Store update after doing L and R - retuns false if no room
2107bool CoinAbcTypeFactorization::storeFT(
2108#if ABC_SMALL < 3
2109  const
2110#endif
2111  CoinIndexedVector *updateFT)
2112{
2113#if ABC_SMALL < 3
2114  CoinSimplexInt numberNonZero = updateFT->getNumElements();
2115#else
2116  CoinSimplexInt numberNonZero = numberRows_;
2117#endif
2118#ifndef ABC_USE_FUNCTION_POINTERS
2119  if (lengthAreaU_ >= lastEntryByColumnU_ + numberNonZero) {
2120#else
2121  if (lengthAreaUPlus_ >= lastEntryByColumnUPlus_ + (3 * numberNonZero + 1) / 2) {
2122    scatterStruct &scatter = scatterUColumn()[numberRows_];
2123    CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
2124#endif
2125#if ABC_SMALL < 3
2126    const CoinFactorizationDouble *COIN_RESTRICT region = denseVector(updateFT);
2127    const CoinSimplexInt *COIN_RESTRICT regionIndex = updateFT->getIndices();
2128#else
2129    CoinSimplexDouble *COIN_RESTRICT region = updateFT->denseVector();
2130#endif
2131#ifndef ABC_USE_FUNCTION_POINTERS
2132    CoinBigIndex *COIN_RESTRICT startColumnU = startColumnUAddress_;
2133    //we are going to save at end of U
2134    startColumnU[numberRows_] = lastEntryByColumnU_;
2135    CoinSimplexInt *COIN_RESTRICT putIndex = indexRowUAddress_ + lastEntryByColumnU_;
2136    CoinFactorizationDouble *COIN_RESTRICT putElement = elementUAddress_ + lastEntryByColumnU_;
2137#else
2138    scatter.offset = lastEntryByColumnUPlus_;
2139    CoinFactorizationDouble *COIN_RESTRICT putElement = elementUColumnPlus + lastEntryByColumnUPlus_;
2140    CoinSimplexInt *COIN_RESTRICT putIndex = reinterpret_cast< CoinSimplexInt * >(putElement + numberNonZero);
2141#endif
2142#if ABC_SMALL < 3
2143    for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
2144      CoinSimplexInt indexValue = regionIndex[i];
2145      CoinSimplexDouble value = region[indexValue];
2146      putElement[i] = value;
2147      putIndex[i] = indexValue;
2148    }
2149#else
2150    numberNonZero = 0;
2151    for (CoinSimplexInt i = 0; i < numberRows_; i++) {
2152      CoinExponent expValue = ABC_EXPONENT(region[i]);
2153      if (expValue) {
2154        if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
2155          putElement[numberNonZero] = region[i];
2156          putIndex[numberNonZero++] = i;
2157        } else {
2158          region[i] = 0.0;
2159        }
2160      }
2161    }
2162#endif
2163#ifndef ABC_USE_FUNCTION_POINTERS
2164    numberInColumnAddress_[numberRows_] = numberNonZero;
2165    lastEntryByColumnU_ += numberNonZero;
2166#else
2167    scatter.number = numberNonZero;
2168#endif
2169    return true;
2170  } else {
2171    return false;
2172  }
2173}
2174CoinSimplexInt CoinAbcTypeFactorization::updateColumnFT(CoinIndexedVector &regionSparseX)
2175{
2176  CoinIndexedVector *regionSparse = &regionSparseX;
2177  CoinSimplexInt numberNonZero = regionSparse->getNumElements();
2178  toLongArray(regionSparse, 0);
2179#ifdef ABC_ORDERED_FACTORIZATION
2180  // Permute in for Ftran
2181  permuteInForFtran(regionSparseX);
2182#endif
2183  //  ******* L
2184  //printf("a\n");
2185  //regionSparse->checkClean();
2186  updateColumnL(regionSparse
2187#if ABC_SMALL < 2
2188    ,
2189    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2190#endif
2191  );
2192  //printf("b\n");
2193  //regionSparse->checkClean();
2194  //row bits here
2195  updateColumnR(regionSparse
2196#if ABC_SMALL < 2
2197    ,
2198    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2199#endif
2200  );
2201
2202  //printf("c\n");
2203  //regionSparse->checkClean();
2204  bool doFT = storeFT(regionSparse);
2205  //printf("d\n");
2206  //regionSparse->checkClean();
2207  //update counts
2208  //  ******* U
2209  updateColumnU(regionSparse
2210#if ABC_SMALL < 2
2211    ,
2212    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2213#endif
2214  );
2215  //printf("e\n");
2216#if ABC_DEBUG
2217  regionSparse->checkClean();
2218#endif
2219  numberNonZero = regionSparse->getNumElements();
2220  // will be negative if no room
2221  if (doFT)
2222    return numberNonZero;
2223  else
2224    return -numberNonZero;
2225}
2226void CoinAbcTypeFactorization::updateColumnFT(CoinIndexedVector &regionSparseX,
2227  CoinIndexedVector &partialUpdate,
2228  int whichSparse)
2229{
2230  CoinIndexedVector *regionSparse = &regionSparseX;
2231  CoinSimplexInt numberNonZero = regionSparse->getNumElements();
2232  toLongArray(regionSparse, whichSparse);
2233#ifdef ABC_ORDERED_FACTORIZATION
2234  // Permute in for Ftran
2235  permuteInForFtran(regionSparseX);
2236#endif
2237  //  ******* L
2238  //printf("a\n");
2239  //regionSparse->checkClean();
2240  updateColumnL(regionSparse
2241#if ABC_SMALL < 2
2242    ,
2243    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2244#endif
2245#if ABC_PARALLEL
2246      ,
2247    whichSparse
2248#endif
2249  );
2250  //printf("b\n");
2251  //regionSparse->checkClean();
2252  //row bits here
2253  updateColumnR(regionSparse
2254#if ABC_SMALL < 2
2255    ,
2256    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2257#endif
2258#if ABC_PARALLEL
2259      ,
2260    whichSparse
2261#endif
2262  );
2263  numberNonZero = regionSparse->getNumElements();
2264  CoinSimplexInt *COIN_RESTRICT indexSave = partialUpdate.getIndices();
2265  CoinSimplexDouble *COIN_RESTRICT regionSave = partialUpdate.denseVector();
2266  partialUpdate.setNumElements(numberNonZero);
2267  memcpy(indexSave, regionSparse->getIndices(), numberNonZero * sizeof(CoinSimplexInt));
2268  partialUpdate.setPackedMode(false);
2269  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
2270  for (int i = 0; i < numberNonZero; i++) {
2271    int iRow = indexSave[i];
2272    regionSave[iRow] = region[iRow];
2273  }
2274  //  ******* U
2275  updateColumnU(regionSparse
2276#if ABC_SMALL < 2
2277    ,
2278    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2279#endif
2280#if ABC_PARALLEL
2281      ,
2282    whichSparse
2283#endif
2284  );
2285  //printf("e\n");
2286#if ABC_DEBUG
2287  regionSparse->checkClean();
2288#endif
2289}
2290int CoinAbcTypeFactorization::updateColumnFTPart1(CoinIndexedVector &regionSparse)
2291{
2292  toLongArray(&regionSparse, 0);
2293#ifdef ABC_ORDERED_FACTORIZATION
2294  // Permute in for Ftran
2295  permuteInForFtran(regionSparse);
2296#endif
2297  //CoinSimplexInt numberNonZero=regionSparse.getNumElements();
2298  //  ******* L
2299  //regionSparse.checkClean();
2300  updateColumnL(&regionSparse
2301#if ABC_SMALL < 2
2302    ,
2303    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2304#endif
2305#if ABC_PARALLEL
2306      ,
2307    2
2308#endif
2309  );
2310  //regionSparse.checkClean();
2311  //row bits here
2312  updateColumnR(&regionSparse
2313#if ABC_SMALL < 2
2314    ,
2315    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2316#endif
2317#if ABC_PARALLEL
2318      ,
2319    2
2320#endif
2321  );
2322
2323  //regionSparse.checkClean();
2324  bool doFT = storeFT(&regionSparse);
2325  // will be negative if no room
2326  if (doFT)
2327    return 1;
2328  else
2329    return -1;
2330}
2331void CoinAbcTypeFactorization::updateColumnFTPart2(CoinIndexedVector &regionSparse)
2332{
2333  toLongArray(&regionSparse, 0);
2334  //CoinSimplexInt numberNonZero=regionSparse.getNumElements();
2335  //  ******* U
2336  updateColumnU(&regionSparse
2337#if ABC_SMALL < 2
2338    ,
2339    reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2340#endif
2341#if ABC_PARALLEL
2342      ,
2343    2
2344#endif
2345  );
2346#if ABC_DEBUG
2347  regionSparse.checkClean();
2348#endif
2349}
2350/* Updates one column (FTRAN) */
2351void CoinAbcTypeFactorization::updateColumnCpu(CoinIndexedVector &regionSparse, int whichCpu) const
2352{
2353  toLongArray(&regionSparse, whichCpu);
2354#ifdef ABC_ORDERED_FACTORIZATION
2355  // Permute in for Ftran
2356  permuteInForFtran(regionSparse);
2357#endif
2358  //  ******* L
2359  updateColumnL(&regionSparse
2360#if ABC_SMALL < 2
2361    ,
2362    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
2363#endif
2364#if ABC_PARALLEL
2365      ,
2366    whichCpu
2367#endif
2368  );
2369  //row bits here
2370  updateColumnR(&regionSparse
2371#if ABC_SMALL < 2
2372    ,
2373    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
2374#endif
2375#if ABC_PARALLEL
2376      ,
2377    whichCpu
2378#endif
2379  );
2380  //update counts
2381  //  ******* U
2382  updateColumnU(&regionSparse
2383#if ABC_SMALL < 2
2384    ,
2385    reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
2386#endif
2387#if ABC_PARALLEL
2388      ,
2389    whichCpu
2390#endif
2391  );
2392  fromLongArray(whichCpu);
2393#if ABC_DEBUG
2394  regionSparse.checkClean();
2395#endif
2396}
2397/* Updates one column (BTRAN) */
2398void CoinAbcTypeFactorization::updateColumnTransposeCpu(CoinIndexedVector &regionSparse, int whichCpu) const
2399{
2400  toLongArray(&regionSparse, whichCpu);
2401
2402  CoinSimplexDouble *COIN_RESTRICT region = regionSparse.denseVector();
2403  CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse.getIndices();
2404  CoinSimplexInt numberNonZero = regionSparse.getNumElements();
2405  //#if COIN_BIG_DOUBLE==1
2406  //const CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegion_.array()+1;
2407  //#else
2408  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
2409  //#endif
2410
2411#ifndef ABC_ORDERED_FACTORIZATION
2412  // can I move this
2413#ifndef INLINE_IT3
2414  for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
2415    CoinSimplexInt iRow = regionIndex[i];
2416    CoinSimplexDouble value = region[iRow];
2417    value *= pivotRegion[iRow];
2418    region[iRow] = value;
2419  }
2420#else
2421  multiplyIndexed(numberNonZero, pivotRegion,
2422    regionIndex, region);
2423#endif
2424#else
2425  // Permute in for Btran
2426  permuteInForBtranAndMultiply(regionSparse);
2427#endif
2428  //  ******* U
2429  // Apply pivot region - could be combined for speed
2430  // Can only combine/move inside vector for sparse
2431  CoinSimplexInt smallestIndex = pivotLinkedForwardsAddress_[-1];
2432#if ABC_SMALL < 2
2433  // copy of code inside transposeU
2434  bool goSparse = false;
2435#else
2436#define goSparse false
2437#endif
2438#if ABC_SMALL < 2
2439  // Guess at number at end
2440  if (gotUCopy()) {
2441    assert(btranAverageAfterU_);
2442    CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * btranAverageAfterU_ * twiddleBtranFactor1());
2443    if (newNumber < sparseThreshold_)
2444      goSparse = true;
2445  }
2446#endif
2447  if (numberNonZero < 40 && (numberNonZero << 4) < numberRows_ && !goSparse) {
2448    CoinSimplexInt *COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
2449    CoinSimplexInt smallest = numberRowsExtra_;
2450    for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
2451      CoinSimplexInt iRow = regionIndex[j];
2452      if (pivotRowForward[iRow] < smallest) {
2453        smallest = pivotRowForward[iRow];
2454        smallestIndex = iRow;
2455      }
2456    }
2457  }
2458  updateColumnTransposeU(&regionSparse, smallestIndex
2459#if ABC_SMALL < 2
2460    ,
2461    reinterpret_cast< CoinAbcStatistics & >(btranCountInput_)
2462#endif
2463#if ABC_PARALLEL
2464      ,
2465    whichCpu
2466#endif
2467  );
2468  //row bits here
2469  updateColumnTransposeR(&regionSparse
2470#if ABC_SMALL < 2
2471    ,
2472    reinterpret_cast< CoinAbcStatistics & >(btranCountInput_)
2473#endif
2474  );
2475  //  ******* L
2476  updateColumnTransposeL(&regionSparse
2477#if ABC_SMALL < 2
2478    ,
2479    reinterpret_cast< CoinAbcStatistics & >(btranCountInput_)
2480#endif
2481#if ABC_PARALLEL
2482      ,
2483    whichCpu
2484#endif
2485  );
2486#if ABC_SMALL < 3
2487#ifdef ABC_ORDERED_FACTORIZATION
2488  // Permute out for Btran
2489  permuteOutForBtran(regionSparse);
2490#endif
2491#if ABC_DEBUG
2492  regionSparse.checkClean();
2493#endif
2494  numberNonZero = regionSparse.getNumElements();
2495#else
2496  numberNonZero = 0;
2497  for (CoinSimplexInt i = 0; i < numberRows_; i++) {
2498    CoinExponent expValue = ABC_EXPONENT(region[i]);
2499    if (expValue) {
2500      if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
2501        regionIndex[numberNonZero++] = i;
2502      } else {
2503        region[i] = 0.0;
2504      }
2505    }
2506  }
2507  regionSparse.setNumElements(numberNonZero);
2508#endif
2509}
2510#if ABC_SMALL < 2
2511//  getRowSpaceIterate.  Gets space for one Row with given length
2512//may have to do compression  (returns true)
2513//also moves existing vector
2514bool CoinAbcTypeFactorization::getRowSpaceIterate(CoinSimplexInt iRow,
2515  CoinSimplexInt extraNeeded)
2516{
2517  const CoinSimplexInt *COIN_RESTRICT numberInRow = numberInRowAddress_;
2518  CoinSimplexInt number = numberInRow[iRow];
2519  CoinBigIndex *COIN_RESTRICT startRow = startRowUAddress_;
2520  CoinSimplexInt *COIN_RESTRICT indexColumnU = indexColumnUAddress_;
2521#if CONVERTROW
2522  CoinBigIndex *COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
2523#if CONVERTROW > 2
2524  CoinBigIndex *COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
2525#endif
2526#endif
2527  CoinFactorizationDouble *COIN_RESTRICT elementRowU = elementRowUAddress_;
2528  CoinBigIndex space = lengthAreaU_ - lastEntryByRowU_;
2529  CoinSimplexInt *COIN_RESTRICT nextRow = nextRowAddress_;
2530  CoinSimplexInt *COIN_RESTRICT lastRow = lastRowAddress_;
2531  if (space < extraNeeded + number + 2) {
2532    //compression
2533    CoinSimplexInt iRow = nextRow[numberRows_];
2534    CoinBigIndex put = 0;
2535    while (iRow != numberRows_) {
2536      //move
2537      CoinBigIndex get = startRow[iRow];
2538      CoinBigIndex getEnd = startRow[iRow] + numberInRow[iRow];
2539
2540      startRow[iRow] = put;
2541      CoinBigIndex i;
2542      for (i = get; i < getEnd; i++) {
2543        indexColumnU[put] = indexColumnU[i];
2544#if CONVERTROW
2545        convertRowToColumn[put] = convertRowToColumn[i];
2546#if CONVERTROW > 2
2547        convertColumnToRow[i] = convertColumnToRow[put];
2548#endif
2549#endif
2550        elementRowU[put] = elementRowU[i];
2551        put++;
2552      }
2553      iRow = nextRow[iRow];
2554    } /* endwhile */
2555    numberCompressions_++;
2556    lastEntryByRowU_ = put;
2557    space = lengthAreaU_ - put;
2558    if (space < extraNeeded + number + 2) {
2559      //need more space
2560      //if we can allocate bigger then do so and copy
2561      //if not then return so code can start again
2562      status_ = -99;
2563      return false;
2564    }
2565  }
2566  CoinBigIndex put = lastEntryByRowU_;
2567  CoinSimplexInt next = nextRow[iRow];
2568  CoinSimplexInt last = lastRow[iRow];
2569
2570  //out
2571  nextRow[last] = next;
2572  lastRow[next] = last;
2573  //in at end
2574  last = lastRow[numberRows_];
2575  nextRow[last] = iRow;
2576  lastRow[numberRows_] = iRow;
2577  lastRow[iRow] = last;
2578  nextRow[iRow] = numberRows_;
2579  //move
2580  CoinBigIndex get = startRow[iRow];
2581  startRow[iRow] = put;
2582  while (number) {
2583    number--;
2584    indexColumnU[put] = indexColumnU[get];
2585#if CONVERTROW
2586    convertRowToColumn[put] = convertRowToColumn[get];
2587#if CONVERTROW > 2
2588    convertColumnToRow[get] = convertColumnToRow[put];
2589#endif
2590#endif
2591    elementRowU[put] = elementRowU[get];
2592    put++;
2593    get++;
2594  } /* endwhile */
2595  //add four for luck
2596  lastEntryByRowU_ = put + extraNeeded + 4;
2597  return true;
2598}
2599#endif
2600void CoinAbcTypeFactorization::printRegion(const CoinIndexedVector &vector, const char *where) const
2601{
2602  //return;
2603  CoinSimplexInt numberNonZero = vector.getNumElements();
2604  //CoinSimplexInt * COIN_RESTRICT regionIndex = vector.getIndices (  );
2605  const CoinSimplexDouble *COIN_RESTRICT region = vector.denseVector();
2606  int n = 0;
2607  for (int i = 0; i < numberRows_; i++) {
2608    if (region[i])
2609      n++;
2610  }
2611  printf("==== %d nonzeros (%d in count) %s ====\n", n, numberNonZero, where);
2612  for (int i = 0; i < numberRows_; i++) {
2613    if (region[i])
2614      printf("%d %g\n", i, region[i]);
2615  }
2616  printf("====            %s ====\n", where);
2617}
2618void CoinAbcTypeFactorization::unpack(CoinIndexedVector *regionFrom,
2619  CoinIndexedVector *regionTo) const
2620{
2621  // move
2622  CoinSimplexInt *COIN_RESTRICT regionIndex = regionTo->getIndices();
2623  CoinSimplexInt numberNonZero = regionFrom->getNumElements();
2624  CoinSimplexInt *COIN_RESTRICT index = regionFrom->getIndices();
2625  CoinSimplexDouble *COIN_RESTRICT region = regionTo->denseVector();
2626  CoinSimplexDouble *COIN_RESTRICT array = regionFrom->denseVector();
2627  for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
2628    CoinSimplexInt iRow = index[j];
2629    CoinSimplexDouble value = array[j];
2630    array[j] = 0.0;
2631    region[iRow] = value;
2632    regionIndex[j] = iRow;
2633  }
2634  regionTo->setNumElements(numberNonZero);
2635}
2636void CoinAbcTypeFactorization::pack(CoinIndexedVector *regionFrom,
2637  CoinIndexedVector *regionTo) const
2638{
2639  // move
2640  CoinSimplexInt *COIN_RESTRICT regionIndex = regionFrom->getIndices();
2641  CoinSimplexInt numberNonZero = regionFrom->getNumElements();
2642  CoinSimplexInt *COIN_RESTRICT index = regionTo->getIndices();
2643  CoinSimplexDouble *COIN_RESTRICT region = regionFrom->denseVector();
2644  CoinSimplexDouble *COIN_RESTRICT array = regionTo->denseVector();
2645  for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
2646    CoinSimplexInt iRow = regionIndex[j];
2647    CoinSimplexDouble value = region[iRow];
2648    region[iRow] = 0.0;
2649    array[j] = value;
2650    index[j] = iRow;
2651  }
2652  regionTo->setNumElements(numberNonZero);
2653}
2654// Set up addresses from arrays
2655void CoinAbcTypeFactorization::doAddresses()
2656{
2657  pivotColumnAddress_ = pivotColumn_.array();
2658  permuteAddress_ = permute_.array();
2659  pivotRegionAddress_ = pivotRegion_.array() + 1;
2660  elementUAddress_ = elementU_.array();
2661  indexRowUAddress_ = indexRowU_.array();
2662  numberInColumnAddress_ = numberInColumn_.array();
2663  numberInColumnPlusAddress_ = numberInColumnPlus_.array();
2664#ifdef ABC_USE_FUNCTION_POINTERS
2665  scatterPointersUColumnAddress_ = reinterpret_cast< scatterStruct * >(scatterUColumn_.array());
2666#endif
2667  startColumnUAddress_ = startColumnU_.array();
2668#if CONVERTROW
2669  convertRowToColumnUAddress_ = convertRowToColumnU_.array();
2670#if CONVERTROW > 1
2671  convertColumnToRowUAddress_ = convertColumnToRowU_.array();
2672#endif
2673#endif
2674#if ABC_SMALL < 2
2675  elementRowUAddress_ = elementRowU_.array();
2676#endif
2677  startRowUAddress_ = startRowU_.array();
2678  numberInRowAddress_ = numberInRow_.array();
2679  indexColumnUAddress_ = indexColumnU_.array();
2680  firstCountAddress_ = firstCount_.array();
2681  nextCountAddress_ = nextCount();
2682  lastCountAddress_ = lastCount();
2683  nextColumnAddress_ = nextColumn_.array();
2684  lastColumnAddress_ = lastColumn_.array();
2685  nextRowAddress_ = nextRow_.array();
2686  lastRowAddress_ = lastRow_.array();
2687  saveColumnAddress_ = saveColumn_.array();
2688  //saveColumnAddress2_ = saveColumn_.array()+numberRows_;
2689  elementLAddress_ = elementL_.array();
2690  indexRowLAddress_ = indexRowL_.array();
2691  startColumnLAddress_ = startColumnL_.array();
2692#if ABC_SMALL < 2
2693  startRowLAddress_ = startRowL_.array();
2694#endif
2695  pivotLinkedBackwardsAddress_ = pivotLinkedBackwards();
2696  pivotLinkedForwardsAddress_ = pivotLinkedForwards();
2697  pivotLOrderAddress_ = pivotLOrder();
2698  startColumnRAddress_ = startColumnR();
2699  // next two are recomputed cleanup
2700  elementRAddress_ = elementL_.array() + lengthL_;
2701  indexRowRAddress_ = indexRowL_.array() + lengthL_;
2702#if ABC_SMALL < 2
2703  indexColumnLAddress_ = indexColumnL_.array();
2704  elementByRowLAddress_ = elementByRowL_.array();
2705#endif
2706#if ABC_DENSE_CODE
2707  denseAreaAddress_ = denseArea_.array();
2708#endif
2709  workAreaAddress_ = workArea_.array();
2710  workArea2Address_ = workArea2_.array();
2711#if ABC_SMALL < 2
2712  sparseAddress_ = sparse_.array();
2713  CoinAbcStack *stackList = reinterpret_cast< CoinAbcStack * >(sparseAddress_);
2714  listAddress_ = reinterpret_cast< CoinSimplexInt * >(stackList + numberRows_);
2715  markRowAddress_ = reinterpret_cast< CoinCheckZero * >(listAddress_ + numberRows_);
2716#endif
2717}
2718#endif
2719
2720/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2721*/
Note: See TracBrowser for help on using the repository browser.