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

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