source: trunk/Clp/src/AbcSimplexParallel.cpp @ 1878

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

File size: 96.4 KB
Line 
1/* $Id: AbcSimplexTemp.cpp 1785 2011-08-25 10:17:49Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5#include "CoinPragma.hpp"
6
7#include <math.h>
8#if CLP_HAS_ABC
9#include "CoinHelperFunctions.hpp"
10#include "CoinAbcHelperFunctions.hpp"
11#include "AbcSimplexDual.hpp"
12#include "ClpEventHandler.hpp"
13#include "AbcSimplexFactorization.hpp"
14#include "CoinPackedMatrix.hpp"
15#include "CoinIndexedVector.hpp"
16#include "CoinFloatEqual.hpp"
17#include "AbcDualRowDantzig.hpp"
18#include "ClpMessage.hpp"
19#include "ClpLinearObjective.hpp"
20#include <cfloat>
21#include <cassert>
22#include <string>
23#include <stdio.h>
24#include <iostream>
25//#undef AVX2
26#define _mm256_broadcast_sd(x) static_cast<__m256d> (__builtin_ia32_vbroadcastsd256 (x))
27#define _mm256_load_pd(x) *(__m256d *)(x)
28#define _mm256_store_pd (s, x)  *((__m256d *)s)=x
29//#define ABC_DEBUG 2
30/* Reasons to come out:
31   -1 iterations etc
32   -2 inaccuracy
33   -3 slight inaccuracy (and done iterations)
34   +0 looks optimal (might be unbounded - but we will investigate)
35   +1 looks infeasible
36   +3 max iterations
37*/
38int
39AbcSimplexDual::whileIteratingSerial()
40{
41  /* arrays
42     0 - to get tableau row and then for weights update
43     1 - tableau column
44     2 - for flip
45     3 - actual tableau row
46  */
47#ifdef CLP_DEBUG
48  int debugIteration = -1;
49#endif
50  // if can't trust much and long way from optimal then relax
51  //if (largestPrimalError_ > 10.0)
52  //abcFactorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
53  //else
54  //abcFactorization_->relaxAccuracyCheck(1.0);
55  // status stays at -1 while iterating, >=0 finished, -2 to invert
56  // status -3 to go to top without an invert
57  int returnCode = -1;
58#define DELAYED_UPDATE
59  arrayForBtran_=0;
60  setUsedArray(arrayForBtran_);
61  arrayForFtran_=1;
62  setUsedArray(arrayForFtran_);
63  arrayForFlipBounds_=2;
64  setUsedArray(arrayForFlipBounds_);
65  arrayForTableauRow_=3;
66  setUsedArray(arrayForTableauRow_);
67  arrayForDualColumn_=4;
68  setUsedArray(arrayForDualColumn_);
69  arrayForReplaceColumn_=4; //4;
70  arrayForFlipRhs_=5;
71  setUsedArray(arrayForFlipRhs_);
72  dualPivotRow();
73  lastPivotRow_=pivotRow_;
74  if (pivotRow_ >= 0) {
75    // we found a pivot row
76    createDualPricingVectorSerial();
77    do {
78#if ABC_DEBUG
79      checkArrays();
80#endif
81      /*
82        -1 - just move dual in values pass
83        0 - take iteration
84        1 - don't take but continue
85        2 - don't take and break
86      */
87      stateOfIteration_=0;
88      returnCode=-1;
89      // put row of tableau in usefulArray[arrayForTableauRow_]
90      /*
91        Could
92        a) copy [arrayForBtran_] and start off updateWeights earlier
93        b) break transposeTimes into two and do after slack part
94        c) also think if cleaner to go all way with update but just don't do final part
95      */
96      //upperTheta=COIN_DBL_MAX;
97      double saveAcceptable=acceptablePivot_;
98      if (largestPrimalError_>1.0e-5||largestDualError_>1.0e-5) {
99        if (!abcFactorization_->pivots())
100          acceptablePivot_*=1.0e2;
101        else if (abcFactorization_->pivots()<5)
102          acceptablePivot_*=1.0e1;
103      }
104      dualColumn1();
105      acceptablePivot_=saveAcceptable;
106      if ((stateOfProblem_&VALUES_PASS)!=0) {
107        // see if can just move dual
108        if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
109          stateOfIteration_=-1;
110        }
111      }
112      //usefulArray_[arrayForTableauRow_].sortPacked();
113      //usefulArray_[arrayForTableauRow_].print();
114      //usefulArray_[arrayForDualColumn_].setPackedMode(true);
115      //usefulArray_[arrayForDualColumn].sortPacked();
116      //usefulArray_[arrayForDualColumn].print();
117      if (!stateOfIteration_) {
118        // get sequenceIn_
119        dualPivotColumn();
120        if (sequenceIn_ >= 0) {
121          // normal iteration
122          // update the incoming column (and do weights)
123          getTableauColumnFlipAndStartReplaceSerial();
124        } else if (stateOfIteration_!=-1) {
125          stateOfIteration_=2;
126        }
127      }
128      if (!stateOfIteration_) {
129        //      assert (stateOfIteration_!=1);
130        int whatNext=0;
131        whatNext = housekeeping();
132        if (whatNext == 1) {
133          problemStatus_ = -2; // refactorize
134        } else if (whatNext == 2) {
135          // maximum iterations or equivalent
136          problemStatus_ = 3;
137          returnCode = 3;
138          stateOfIteration_=2;
139        }
140        if (problemStatus_==-1) { 
141          replaceColumnPart3();
142          //clearArrays(arrayForReplaceColumn_);
143#if ABC_DEBUG
144          checkArrays();
145#endif
146            updateDualsInDual();
147            abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
148                                                             movement_);
149#if ABC_DEBUG
150          checkArrays();
151#endif
152        } else {
153          abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
154          //clearArrays(arrayForBtran_);
155          //clearArrays(arrayForFtran_);
156        }
157      } else {
158        if (stateOfIteration_<0) {
159          // move dual in dual values pass
160          theta_=abcDj_[sequenceOut_];
161          updateDualsInDual();
162          abcDj_[sequenceOut_]=0.0;
163          sequenceOut_=-1;
164        }
165        // clear all arrays
166        clearArrays(-1);
167        //sequenceIn_=-1;
168        //sequenceOut_=-1;
169      }
170      // Check event
171      {
172        int status = eventHandler_->event(ClpEventHandler::endOfIteration);
173        if (status >= 0) {
174          problemStatus_ = 5;
175          secondaryStatus_ = ClpEventHandler::endOfIteration;
176          returnCode = 4;
177          stateOfIteration_=2;
178        }
179      }
180      // at this stage sequenceIn_ may be <0
181      if (sequenceIn_<0&&sequenceOut_>=0) {
182        // no incoming column is valid
183        returnCode=noPivotColumn();
184      }
185      if (stateOfIteration_==2) {
186        sequenceIn_=-1;
187        break;
188      }
189      swapPrimalStuff();
190      if (problemStatus_!=-1) {
191        break;
192      }
193      // dualRow will go to virtual row pivot choice algorithm
194      // make sure values pass off if it should be
195      // use Btran array and clear inside dualPivotRow (if used)
196      int lastSequenceOut=sequenceOut_;
197      int lastDirectionOut=directionOut_;
198      dualPivotRow();
199      lastPivotRow_=pivotRow_;
200      if (pivotRow_ >= 0) {
201        // we found a pivot row
202        // create as packed
203        createDualPricingVectorSerial();
204        swapDualStuff(lastSequenceOut,lastDirectionOut);
205        // next pivot
206      } else {
207        // no pivot row
208        clearArrays(-1);
209        returnCode=noPivotRow();
210        break;
211      }
212    } while (problemStatus_ == -1);
213    usefulArray_[arrayForTableauRow_].compact();
214#if ABC_DEBUG
215    checkArrays();
216#endif
217  } else {
218    // no pivot row on first try
219    clearArrays(-1);
220    returnCode=noPivotRow();
221  }
222  //printStuff();
223  clearArrays(-1);
224  //if (pivotRow_==-100)
225  //returnCode=-100; // end of values pass
226  return returnCode;
227}
228// Create dual pricing vector
229void 
230AbcSimplexDual::createDualPricingVectorSerial()
231{
232#ifndef NDEBUG
233#if ABC_NORMAL_DEBUG>3
234    checkArrays();
235#endif
236#endif
237  // we found a pivot row
238  if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
239    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
240      << pivotRow_
241      << CoinMessageEol;
242  }
243  // check accuracy of weights
244  abcDualRowPivot_->checkAccuracy();
245  // get sign for finding row of tableau
246  // create as packed
247  usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
248  abcFactorization_->updateColumnTranspose(usefulArray_[arrayForBtran_]);
249  sequenceIn_ = -1;
250}
251void
252AbcSimplexDual::getTableauColumnPart1Serial()
253{
254#if ABC_DEBUG
255  {
256    const double * work=usefulArray_[arrayForTableauRow_].denseVector();
257    int number=usefulArray_[arrayForTableauRow_].getNumElements();
258    const int * which=usefulArray_[arrayForTableauRow_].getIndices();
259    for (int i=0;i<number;i++) {
260      if (which[i]==sequenceIn_) {
261        assert (alpha_==work[i]);
262        break;
263      }
264    }
265  }
266#endif
267  //int returnCode=0;
268  // update the incoming column
269  unpack(usefulArray_[arrayForFtran_]);
270  // Array[0] may be needed until updateWeights2
271  // and update dual weights (can do in parallel - with extra array)
272  alpha_ = abcDualRowPivot_->updateWeights1(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
273}
274int 
275AbcSimplexDual::getTableauColumnFlipAndStartReplaceSerial()
276{
277  // move checking stuff down into called functions
278  int numberFlipped;
279  numberFlipped=flipBounds();
280  // update the incoming column
281  getTableauColumnPart1Serial();
282  checkReplacePart1();
283  //checkReplacePart1a();
284  //checkReplacePart1b();
285  getTableauColumnPart2();
286  //usefulArray_[arrayForTableauRow_].compact();
287  // returns 3 if skip this iteration and re-factorize
288/*
289  return code
290  0 - OK
291  2 - flag something and skip
292  3 - break and re-factorize
293  4 - break and say infeasible
294 */
295  int returnCode=0;
296  // amount primal will move
297  movement_ = -dualOut_ * directionOut_ / alpha_;
298  // see if update stable
299#if ABC_NORMAL_DEBUG>3
300  if ((handler_->logLevel() & 32)) {
301    double ft=ftAlpha_*abcFactorization_->pivotRegion()[pivotRow_];
302    double diff1=fabs(alpha_-btranAlpha_);
303    double diff2=fabs(fabs(alpha_)-fabs(ft));
304    double diff3=fabs(fabs(ft)-fabs(btranAlpha_));
305    double largest=CoinMax(CoinMax(diff1,diff2),diff3);
306    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n", 
307           btranAlpha_, alpha_,ft,largest);
308    if (largest>0.001*fabs(alpha_)) {
309      printf("bad\n");
310    }
311  }
312#endif
313  int numberPivots=abcFactorization_->pivots();
314  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
315  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
316  // if can't trust much and long way from optimal then relax
317  if (largestPrimalError_ > 10.0)
318    checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
319  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
320      fabs(btranAlpha_ - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
321    handler_->message(CLP_DUAL_CHECK, messages_)
322      << btranAlpha_
323      << alpha_
324      << CoinMessageEol;
325    // see with more relaxed criterion
326    double test;
327    if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
328      test = 1.0e-1 * fabs(alpha_);
329    else
330      test = 10.0*checkValue;//1.0e-4 * (1.0 + fabs(alpha_));
331    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
332                     fabs(btranAlpha_ - alpha_) > test);
333    double derror = CoinMin(fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),1.0)*0.9999e7;
334    int error=0;
335    while (derror>1.0) {
336      error++;
337      derror *= 0.1;
338    }
339    int frequency[8]={99999999,100,10,2,1,1,1,1};
340    int newFactorFrequency=CoinMin(forceFactorization_,frequency[error]);
341#if ABC_NORMAL_DEBUG>0
342    if (newFactorFrequency<forceFactorization_)
343      printf("Error of %g after %d pivots old forceFact %d now %d\n",
344             fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),numberPivots,
345             forceFactorization_,newFactorFrequency);
346#endif
347    if (!numberPivots&&fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_))>0.5e-6
348        &&abcFactorization_->pivotTolerance()<0.5)
349      abcFactorization_->saferTolerances (1.0, -1.03);
350    forceFactorization_=newFactorFrequency;
351
352   
353#if ABC_NORMAL_DEBUG>0
354    if (fabs(btranAlpha_ + alpha_) < 1.0e-5*(1.0 + fabs(alpha_))) {
355      printf("diff (%g,%g) %g check %g\n",btranAlpha_,alpha_,fabs(btranAlpha_-alpha_),checkValue*(1.0+fabs(alpha_)));
356    }
357#endif
358    if (numberPivots) {
359      if (needFlag&&numberPivots<10) {
360        // need to reject something
361        assert (sequenceOut_>=0);
362        char x = isColumn(sequenceOut_) ? 'C' : 'R';
363        handler_->message(CLP_SIMPLEX_FLAG, messages_)
364          << x << sequenceWithin(sequenceOut_)
365          << CoinMessageEol;
366#if ABC_NORMAL_DEBUG>0
367        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
368               btranAlpha_, alpha_,numberPivots);
369#endif
370        // Make safer?
371        abcFactorization_->saferTolerances (1.0, -1.03);
372        setFlagged(sequenceOut_);
373        // so can't happen immediately again
374        sequenceOut_=-1;
375        lastBadIteration_ = numberIterations_; // say be more cautious
376      }
377      // Make safer?
378      //if (numberPivots<5)
379      //abcFactorization_->saferTolerances (-0.99, -1.01);
380      problemStatus_ = -2; // factorize now
381      returnCode = -2;
382      stateOfIteration_=2;
383    } else {
384      if (needFlag) {
385        assert (sequenceOut_>=0);
386        // need to reject something
387        char x = isColumn(sequenceOut_) ? 'C' : 'R';
388        handler_->message(CLP_SIMPLEX_FLAG, messages_)
389          << x << sequenceWithin(sequenceOut_)
390          << CoinMessageEol;
391#if ABC_NORMAL_DEBUG>3
392        printf("flag a %g %g\n", btranAlpha_, alpha_);
393#endif
394        setFlagged(sequenceOut_);
395        // so can't happen immediately again
396        sequenceOut_=-1;
397        //abcProgress_.clearBadTimes();
398        lastBadIteration_ = numberIterations_; // say be more cautious
399        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
400          //printf("I think should declare infeasible\n");
401          problemStatus_ = 1;
402          returnCode = 1;
403          stateOfIteration_=2;
404        } else {
405          stateOfIteration_=1;
406        }
407        // Make safer?
408        if (abcFactorization_->pivotTolerance()<0.999&&stateOfIteration_==1) {
409          // change tolerance and re-invert
410          abcFactorization_->saferTolerances (1.0, -1.03);
411          problemStatus_ = -6; // factorize now
412          returnCode = -2;
413          stateOfIteration_=2;
414        }
415      }
416    }
417  }
418  if (!stateOfIteration_) {
419    // check update
420    //int updateStatus =
421/*
422  returns
423  0 - OK
424  1 - take iteration then re-factorize
425  2 - flag something and skip
426  3 - break and re-factorize
427  5 - take iteration then re-factorize because of memory
428 */
429    int status=checkReplace();
430    if (status&&!returnCode)
431      returnCode=status;
432  }
433  //clearArrays(arrayForFlipRhs_);
434  if (stateOfIteration_&&numberFlipped) {
435    //usefulArray_[arrayForTableauRow_].compact();
436    // move variables back
437    flipBack(numberFlipped);
438  }
439  // could choose average of all three
440  return returnCode;
441}
442#if ABC_PARALLEL==1
443/* Reasons to come out:
444   -1 iterations etc
445   -2 inaccuracy
446   -3 slight inaccuracy (and done iterations)
447   +0 looks optimal (might be unbounded - but we will investigate)
448   +1 looks infeasible
449   +3 max iterations
450*/
451int
452AbcSimplexDual::whileIteratingThread()
453{
454  /* arrays
455     0 - to get tableau row and then for weights update
456     1 - tableau column
457     2 - for flip
458     3 - actual tableau row
459  */
460#ifdef CLP_DEBUG
461  int debugIteration = -1;
462#endif
463  // if can't trust much and long way from optimal then relax
464  //if (largestPrimalError_ > 10.0)
465  //abcFactorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
466  //else
467  //abcFactorization_->relaxAccuracyCheck(1.0);
468  // status stays at -1 while iterating, >=0 finished, -2 to invert
469  // status -3 to go to top without an invert
470  int returnCode = -1;
471#define DELAYED_UPDATE
472  arrayForBtran_=0;
473  setUsedArray(arrayForBtran_);
474  arrayForFtran_=1;
475  setUsedArray(arrayForFtran_);
476  arrayForFlipBounds_=2;
477  setUsedArray(arrayForFlipBounds_);
478  arrayForTableauRow_=3;
479  setUsedArray(arrayForTableauRow_);
480  arrayForDualColumn_=4;
481  setUsedArray(arrayForDualColumn_);
482  arrayForReplaceColumn_=4; //4;
483  arrayForFlipRhs_=5;
484  setUsedArray(arrayForFlipRhs_);
485  dualPivotRow();
486  lastPivotRow_=pivotRow_;
487  if (pivotRow_ >= 0) {
488    // we found a pivot row
489    createDualPricingVectorThread();
490    do {
491#if ABC_DEBUG
492      checkArrays();
493#endif
494      /*
495        -1 - just move dual in values pass
496        0 - take iteration
497        1 - don't take but continue
498        2 - don't take and break
499      */
500      stateOfIteration_=0;
501      returnCode=-1;
502      // put row of tableau in usefulArray[arrayForTableauRow_]
503      /*
504        Could
505        a) copy [arrayForBtran_] and start off updateWeights earlier
506        b) break transposeTimes into two and do after slack part
507        c) also think if cleaner to go all way with update but just don't do final part
508      */
509      //upperTheta=COIN_DBL_MAX;
510      double saveAcceptable=acceptablePivot_;
511      if (largestPrimalError_>1.0e-5||largestDualError_>1.0e-5) {
512        if (!abcFactorization_->pivots())
513          acceptablePivot_*=1.0e2;
514        else if (abcFactorization_->pivots()<5)
515          acceptablePivot_*=1.0e1;
516      }
517      dualColumn1();
518      acceptablePivot_=saveAcceptable;
519      if ((stateOfProblem_&VALUES_PASS)!=0) {
520        // see if can just move dual
521        if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
522          stateOfIteration_=-1;
523        }
524      }
525      //usefulArray_[arrayForTableauRow_].sortPacked();
526      //usefulArray_[arrayForTableauRow_].print();
527      //usefulArray_[arrayForDualColumn_].setPackedMode(true);
528      //usefulArray_[arrayForDualColumn].sortPacked();
529      //usefulArray_[arrayForDualColumn].print();
530      if (parallelMode_!=0) {
531        stopStart_=1+32*1; // Just first thread for updateweights
532        startParallelStuff(2);
533      }
534      if (!stateOfIteration_) {
535        // get sequenceIn_
536        dualPivotColumn();
537        if (sequenceIn_ >= 0) {
538          // normal iteration
539          // update the incoming column (and do weights if serial)
540          getTableauColumnFlipAndStartReplaceThread();
541          //usleep(1000);
542        } else if (stateOfIteration_!=-1) {
543          stateOfIteration_=2;
544        }
545      }
546      if (parallelMode_!=0) {
547        // do sync here
548        stopParallelStuff(2);
549      }
550      if (!stateOfIteration_) {
551        //      assert (stateOfIteration_!=1);
552        int whatNext=0;
553        whatNext = housekeeping();
554        if (whatNext == 1) {
555          problemStatus_ = -2; // refactorize
556        } else if (whatNext == 2) {
557          // maximum iterations or equivalent
558          problemStatus_ = 3;
559          returnCode = 3;
560          stateOfIteration_=2;
561        }
562      } else {
563        if (stateOfIteration_<0) {
564          // move dual in dual values pass
565          theta_=abcDj_[sequenceOut_];
566          updateDualsInDual();
567          abcDj_[sequenceOut_]=0.0;
568          sequenceOut_=-1;
569        }
570        // clear all arrays
571        clearArrays(-1);
572      }
573      // at this stage sequenceIn_ may be <0
574      if (sequenceIn_<0&&sequenceOut_>=0) {
575        // no incoming column is valid
576        returnCode=noPivotColumn();
577      }
578      if (stateOfIteration_==2) {
579        sequenceIn_=-1;
580        break;
581      }
582      if (problemStatus_!=-1) {
583        abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
584        swapPrimalStuff();
585        break;
586      }
587#if ABC_DEBUG
588      checkArrays();
589#endif
590      if (stateOfIteration_!=-1) {
591        assert (stateOfIteration_==0); // if 1 why are we here
592        // can do these in parallel
593        if (parallelMode_==0) {
594          replaceColumnPart3();
595          updateDualsInDual();
596          abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
597                                                           movement_);
598        } else {
599          stopStart_=1+32*1;
600          startParallelStuff(3);
601          if (parallelMode_>1) {
602            stopStart_=2+32*2;
603            startParallelStuff(9);
604          } else {
605            replaceColumnPart3();
606          }
607          abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
608                                                           movement_);
609        }
610      }
611#if ABC_DEBUG
612      checkArrays();
613#endif
614      swapPrimalStuff();
615      // dualRow will go to virtual row pivot choice algorithm
616      // make sure values pass off if it should be
617      // use Btran array and clear inside dualPivotRow (if used)
618      int lastSequenceOut=sequenceOut_;
619      int lastDirectionOut=directionOut_;
620      // redo to test on parallelMode_
621      if (parallelMode_>1) {
622        stopStart_=2+32*2; // stop factorization update
623        //stopStart_=3+32*3; // stop all
624        stopParallelStuff(9);
625      }
626      dualPivotRow();
627      lastPivotRow_=pivotRow_;
628      if (pivotRow_ >= 0) {
629        // we found a pivot row
630        // create as packed
631        createDualPricingVectorThread();
632        swapDualStuff(lastSequenceOut,lastDirectionOut);
633        // next pivot
634        // redo to test on parallelMode_
635        if (parallelMode_!=0) {
636          stopStart_=1+32*1; // stop dual update
637          stopParallelStuff(3);
638        }
639      } else {
640        // no pivot row
641        // redo to test on parallelMode_
642        if (parallelMode_!=0) {
643          stopStart_=1+32*1; // stop dual update
644          stopParallelStuff(3);
645        }
646        clearArrays(-1);
647        returnCode=noPivotRow();
648        break;
649      }
650    } while (problemStatus_ == -1);
651    usefulArray_[arrayForTableauRow_].compact();
652#if ABC_DEBUG
653    checkArrays();
654#endif
655  } else {
656    // no pivot row on first try
657    clearArrays(-1);
658    returnCode=noPivotRow();
659  }
660  //printStuff();
661  clearArrays(-1);
662  //if (pivotRow_==-100)
663  //returnCode=-100; // end of values pass
664  return returnCode;
665}
666// Create dual pricing vector
667void 
668AbcSimplexDual::createDualPricingVectorThread()
669{
670#ifndef NDEBUG
671#if ABC_NORMAL_DEBUG>3
672    checkArrays();
673#endif
674#endif
675  // we found a pivot row
676  if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
677    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
678      << pivotRow_
679      << CoinMessageEol;
680  }
681  // check accuracy of weights
682  abcDualRowPivot_->checkAccuracy();
683  // get sign for finding row of tableau
684  // create as packed
685  usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
686  abcFactorization_->updateColumnTranspose(usefulArray_[arrayForBtran_]);
687  sequenceIn_ = -1;
688}
689void
690AbcSimplexDual::getTableauColumnPart1Thread()
691{
692#if ABC_DEBUG
693  {
694    const double * work=usefulArray_[arrayForTableauRow_].denseVector();
695    int number=usefulArray_[arrayForTableauRow_].getNumElements();
696    const int * which=usefulArray_[arrayForTableauRow_].getIndices();
697    for (int i=0;i<number;i++) {
698      if (which[i]==sequenceIn_) {
699        assert (alpha_==work[i]);
700        break;
701      }
702    }
703  }
704#endif
705  //int returnCode=0;
706  // update the incoming column
707  unpack(usefulArray_[arrayForFtran_]);
708  // Array[0] may be needed until updateWeights2
709  // and update dual weights (can do in parallel - with extra array)
710  if (parallelMode_!=0) {
711    abcFactorization_->updateColumnFTPart1(usefulArray_[arrayForFtran_]);
712    // pivot element
713    //alpha_ = usefulArray_[arrayForFtran_].denseVector()[pivotRow_];
714    } else {
715    alpha_ = abcDualRowPivot_->updateWeights1(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
716  }
717}
718int 
719AbcSimplexDual::getTableauColumnFlipAndStartReplaceThread()
720{
721  // move checking stuff down into called functions
722  // threads 2 and 3 are available
723  int numberFlipped;
724
725  if (parallelMode_==0) {
726    numberFlipped=flipBounds();
727    // update the incoming column
728    getTableauColumnPart1Thread();
729    checkReplacePart1();
730    //checkReplacePart1a();
731    //checkReplacePart1b();
732    getTableauColumnPart2();
733  } else {
734    // save stopStart
735    int saveStopStart=stopStart_;
736    if (parallelMode_>1) {
737      // we can flip immediately
738      stopStart_=2+32*2; // just thread 1
739      startParallelStuff(8);
740      // update the incoming column
741      getTableauColumnPart1Thread();
742      stopStart_=4+32*4; // just thread 2
743      startParallelStuff(7);
744      getTableauColumnPart2();
745      stopStart_=6+32*6; // just threads 1 and 2
746      stopParallelStuff(8);
747      //ftAlpha_=threadInfo_[2].result;
748    } else {
749    // just one extra thread - do replace
750    // update the incoming column
751    getTableauColumnPart1Thread();
752    stopStart_=2+32*2; // just thread 1
753    startParallelStuff(7);
754    getTableauColumnPart2();
755    numberFlipped=flipBounds();
756    stopParallelStuff(8);
757    //ftAlpha_=threadInfo_[1].result;
758    }
759    stopStart_=saveStopStart;
760  }
761  //usefulArray_[arrayForTableauRow_].compact();
762  // returns 3 if skip this iteration and re-factorize
763/*
764  return code
765  0 - OK
766  2 - flag something and skip
767  3 - break and re-factorize
768  4 - break and say infeasible
769 */
770  int returnCode=0;
771  // amount primal will move
772  movement_ = -dualOut_ * directionOut_ / alpha_;
773  // see if update stable
774#if ABC_NORMAL_DEBUG>3
775  if ((handler_->logLevel() & 32)) {
776    double ft=ftAlpha_*abcFactorization_->pivotRegion()[pivotRow_];
777    double diff1=fabs(alpha_-btranAlpha_);
778    double diff2=fabs(fabs(alpha_)-fabs(ft));
779    double diff3=fabs(fabs(ft)-fabs(btranAlpha_));
780    double largest=CoinMax(CoinMax(diff1,diff2),diff3);
781    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n", 
782           btranAlpha_, alpha_,ft,largest);
783    if (largest>0.001*fabs(alpha_)) {
784      printf("bad\n");
785    }
786  }
787#endif
788  int numberPivots=abcFactorization_->pivots();
789  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
790  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
791  // if can't trust much and long way from optimal then relax
792  if (largestPrimalError_ > 10.0)
793    checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
794  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
795      fabs(btranAlpha_ - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
796    handler_->message(CLP_DUAL_CHECK, messages_)
797      << btranAlpha_
798      << alpha_
799      << CoinMessageEol;
800    // see with more relaxed criterion
801    double test;
802    if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
803      test = 1.0e-1 * fabs(alpha_);
804    else
805      test = 10.0*checkValue;//1.0e-4 * (1.0 + fabs(alpha_));
806    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
807                     fabs(btranAlpha_ - alpha_) > test);
808    double derror = CoinMin(fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),1.0)*0.9999e7;
809    int error=0;
810    while (derror>1.0) {
811      error++;
812      derror *= 0.1;
813    }
814    int frequency[8]={99999999,100,10,2,1,1,1,1};
815    int newFactorFrequency=CoinMin(forceFactorization_,frequency[error]);
816#if ABC_NORMAL_DEBUG>0
817    if (newFactorFrequency<forceFactorization_)
818      printf("Error of %g after %d pivots old forceFact %d now %d\n",
819             fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),numberPivots,
820             forceFactorization_,newFactorFrequency);
821#endif
822    if (!numberPivots&&fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_))>0.5e-6
823        &&abcFactorization_->pivotTolerance()<0.5)
824      abcFactorization_->saferTolerances (1.0, -1.03);
825    forceFactorization_=newFactorFrequency;
826
827   
828#if ABC_NORMAL_DEBUG>0
829    if (fabs(btranAlpha_ + alpha_) < 1.0e-5*(1.0 + fabs(alpha_))) {
830      printf("diff (%g,%g) %g check %g\n",btranAlpha_,alpha_,fabs(btranAlpha_-alpha_),checkValue*(1.0+fabs(alpha_)));
831    }
832#endif
833    if (numberPivots) {
834      if (needFlag&&numberPivots<10) {
835        // need to reject something
836        assert (sequenceOut_>=0);
837        char x = isColumn(sequenceOut_) ? 'C' : 'R';
838        handler_->message(CLP_SIMPLEX_FLAG, messages_)
839          << x << sequenceWithin(sequenceOut_)
840          << CoinMessageEol;
841#if ABC_NORMAL_DEBUG>0
842        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
843               btranAlpha_, alpha_,numberPivots);
844#endif
845        // Make safer?
846        abcFactorization_->saferTolerances (1.0, -1.03);
847        setFlagged(sequenceOut_);
848        // so can't happen immediately again
849        sequenceOut_=-1;
850        lastBadIteration_ = numberIterations_; // say be more cautious
851      }
852      // Make safer?
853      //if (numberPivots<5)
854      //abcFactorization_->saferTolerances (-0.99, -1.01);
855      problemStatus_ = -2; // factorize now
856      returnCode = -2;
857      stateOfIteration_=2;
858    } else {
859      if (needFlag) {
860        assert (sequenceOut_>=0);
861        // need to reject something
862        char x = isColumn(sequenceOut_) ? 'C' : 'R';
863        handler_->message(CLP_SIMPLEX_FLAG, messages_)
864          << x << sequenceWithin(sequenceOut_)
865          << CoinMessageEol;
866#if ABC_NORMAL_DEBUG>3
867        printf("flag a %g %g\n", btranAlpha_, alpha_);
868#endif
869        setFlagged(sequenceOut_);
870        // so can't happen immediately again
871        sequenceOut_=-1;
872        //abcProgress_.clearBadTimes();
873        lastBadIteration_ = numberIterations_; // say be more cautious
874        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
875          //printf("I think should declare infeasible\n");
876          problemStatus_ = 1;
877          returnCode = 1;
878          stateOfIteration_=2;
879        } else {
880          stateOfIteration_=1;
881        }
882        // Make safer?
883        if (abcFactorization_->pivotTolerance()<0.999&&stateOfIteration_==1) {
884          // change tolerance and re-invert
885          abcFactorization_->saferTolerances (1.0, -1.03);
886          problemStatus_ = -6; // factorize now
887          returnCode = -2;
888          stateOfIteration_=2;
889        }
890      }
891    }
892  }
893  if (!stateOfIteration_) {
894    // check update
895    //int updateStatus =
896/*
897  returns
898  0 - OK
899  1 - take iteration then re-factorize
900  2 - flag something and skip
901  3 - break and re-factorize
902  5 - take iteration then re-factorize because of memory
903 */
904    int status=checkReplace();
905    if (status&&!returnCode)
906      returnCode=status;
907  }
908  //clearArrays(arrayForFlipRhs_);
909  if (stateOfIteration_&&numberFlipped) {
910    //usefulArray_[arrayForTableauRow_].compact();
911    // move variables back
912    flipBack(numberFlipped);
913  }
914  // could choose average of all three
915  return returnCode;
916}
917#endif
918#if ABC_PARALLEL==2
919#if ABC_NORMAL_DEBUG>0
920// for conflicts
921int cilk_conflict=0;
922#endif
923#ifdef EARLY_FACTORIZE
924static int doEarlyFactorization(AbcSimplexDual * dual)
925{
926  int returnCode=cilk_spawn dual->whileIteratingParallel(123456789);
927  CoinIndexedVector & vector = *dual->usefulArray(ABC_NUMBER_USEFUL-1);
928  int status=dual->earlyFactorization()->factorize(dual,vector);
929#if 0
930  // Is this safe
931  if (!status&&true) {
932    // do pivots if there are any
933    int capacity = vector.capacity();
934    capacity--;
935    int numberPivotsStored = vector.getIndices()[capacity];
936    status =
937      dual->earlyFactorization()->replaceColumns(dual,vector,
938                                                    0,numberPivotsStored,false);
939  }
940#endif
941  if (status) {
942    printf("bad early factorization in doEarly - switch off\n");
943    vector.setNumElements(-1);
944  }
945  cilk_sync;
946  return returnCode;
947}
948#endif
949#define MOVE_UPDATE_WEIGHTS
950/* Reasons to come out:
951   -1 iterations etc
952   -2 inaccuracy
953   -3 slight inaccuracy (and done iterations)
954   +0 looks optimal (might be unbounded - but we will investigate)
955   +1 looks infeasible
956   +3 max iterations
957*/
958int
959AbcSimplexDual::whileIteratingCilk()
960{
961  /* arrays
962     0 - to get tableau row and then for weights update
963     1 - tableau column
964     2 - for flip
965     3 - actual tableau row
966  */
967#ifdef CLP_DEBUG
968  int debugIteration = -1;
969#endif
970  // if can't trust much and long way from optimal then relax
971  //if (largestPrimalError_ > 10.0)
972  //abcFactorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
973  //else
974  //abcFactorization_->relaxAccuracyCheck(1.0);
975  // status stays at -1 while iterating, >=0 finished, -2 to invert
976  // status -3 to go to top without an invert
977  int returnCode = -1;
978#define DELAYED_UPDATE
979  arrayForBtran_=0;
980  setUsedArray(arrayForBtran_);
981  arrayForFtran_=1;
982  setUsedArray(arrayForFtran_);
983  arrayForFlipBounds_=2;
984  setUsedArray(arrayForFlipBounds_);
985  arrayForTableauRow_=3;
986  setUsedArray(arrayForTableauRow_);
987  arrayForDualColumn_=4;
988  setUsedArray(arrayForDualColumn_);
989  arrayForReplaceColumn_=6; //4;
990  setUsedArray(arrayForReplaceColumn_);
991#ifndef MOVE_UPDATE_WEIGHTS
992  arrayForFlipRhs_=5;
993  setUsedArray(arrayForFlipRhs_);
994#else
995  arrayForFlipRhs_=0;
996  // use for weights
997  setUsedArray(5);
998#endif
999  dualPivotRow();
1000  lastPivotRow_=pivotRow_;
1001  if (pivotRow_ >= 0) {
1002    // we found a pivot row
1003    createDualPricingVectorCilk();
1004    // ABC_PARALLEL 2
1005#if 0
1006    {
1007      for (int i=0;i<numberRows_;i++) {
1008        int iSequence=abcPivotVariable_[i];
1009        if (lowerBasic_[i]!=lowerSaved_[iSequence]||
1010            upperBasic_[i]!=upperSaved_[iSequence])
1011          printf("basic l %g,%g,%g u %g,%g\n",
1012                 abcLower_[iSequence],lowerSaved_[iSequence],lowerBasic_[i],
1013                 abcUpper_[iSequence],upperSaved_[iSequence],upperBasic_[i]);
1014      }
1015    }
1016#endif
1017#if 0
1018    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1019      int bad=-1;
1020      if (getFakeBound(iSequence)==noFake) {
1021        if (abcLower_[iSequence]!=lowerSaved_[iSequence]||
1022            abcUpper_[iSequence]!=upperSaved_[iSequence])
1023          bad=0;
1024      } else if (getFakeBound(iSequence)==lowerFake) {
1025        if (abcLower_[iSequence]==lowerSaved_[iSequence]||
1026            abcUpper_[iSequence]!=upperSaved_[iSequence])
1027          bad=1;
1028      } else if (getFakeBound(iSequence)==upperFake) {
1029        if (abcLower_[iSequence]!=lowerSaved_[iSequence]||
1030            abcUpper_[iSequence]==upperSaved_[iSequence])
1031          bad=2;
1032      } else {
1033        if (abcLower_[iSequence]==lowerSaved_[iSequence]||
1034            abcUpper_[iSequence]==upperSaved_[iSequence])
1035          bad=3;
1036      }
1037      if (bad>=0)
1038        printf("%d status %d l %g,%g u %g,%g\n",iSequence,internalStatus_[iSequence],
1039               abcLower_[iSequence],lowerSaved_[iSequence],
1040               abcUpper_[iSequence],upperSaved_[iSequence]);
1041    }
1042#endif
1043    int numberPivots = abcFactorization_->maximumPivots();
1044#ifdef EARLY_FACTORIZE
1045    int numberEarly=0;
1046    if (numberPivots>20&&(numberEarly_&0xffff)>5) {
1047      numberEarly=numberEarly_&0xffff;
1048      numberPivots=CoinMax(numberPivots-numberEarly-abcFactorization_->pivots(),numberPivots/2);
1049    }
1050    returnCode=whileIteratingParallel(numberPivots);
1051    if (returnCode==-99) {
1052      if (numberEarly) {
1053        if (!abcEarlyFactorization_)
1054          abcEarlyFactorization_ = new AbcSimplexFactorization(*abcFactorization_);
1055        // create pivot list
1056        CoinIndexedVector & vector = usefulArray_[ABC_NUMBER_USEFUL-1];
1057        vector.checkClear();
1058        int * indices = vector.getIndices();
1059        int capacity = vector.capacity();
1060        int numberNeeded=numberRows_+2*numberEarly+1;
1061        if (capacity<numberNeeded) {
1062          vector.reserve(numberNeeded);
1063          capacity = vector.capacity();
1064        }
1065        int numberBasic=0;
1066        for (int i=0;i<numberRows_;i++) {
1067          int iSequence = abcPivotVariable_[i];
1068          if (iSequence<numberRows_)
1069            indices[numberBasic++]=iSequence;
1070        }
1071        vector.setNumElements(numberRows_-numberBasic);
1072        for (int i=0;i<numberRows_;i++) {
1073          int iSequence = abcPivotVariable_[i];
1074          if (iSequence>=numberRows_)
1075            indices[numberBasic++]=iSequence;
1076        }
1077        assert (numberBasic==numberRows_);
1078        indices[capacity-1]=0;
1079        // could set iterations to -1 for safety
1080#if 0
1081        cilk_spawn doEarlyFactorization(this);
1082        returnCode=whileIteratingParallel(123456789);
1083        cilk_sync;
1084#else
1085        returnCode=doEarlyFactorization(this);
1086#endif
1087      } else {
1088        returnCode=whileIteratingParallel(123456789);
1089      }
1090    }
1091#else
1092    returnCode=whileIteratingParallel(numberPivots);
1093#endif
1094    usefulArray_[arrayForTableauRow_].compact();
1095#if ABC_DEBUG
1096    checkArrays();
1097#endif
1098  } else {
1099    // no pivot row on first try
1100    clearArrays(-1);
1101    returnCode=noPivotRow();
1102  }
1103  //printStuff();
1104  clearArrays(-1);
1105  //if (pivotRow_==-100)
1106  //returnCode=-100; // end of values pass
1107  return returnCode;
1108}
1109// Create dual pricing vector
1110void 
1111AbcSimplexDual::createDualPricingVectorCilk()
1112{
1113#if CILK_CONFLICT>0
1114  if(cilk_conflict&15!=0) {
1115    printf("cilk_conflict %d!\n",cilk_conflict);
1116    cilk_conflict=0;
1117  }
1118#endif
1119#ifndef NDEBUG
1120#if ABC_NORMAL_DEBUG>3
1121    checkArrays();
1122#endif
1123#endif
1124  // we found a pivot row
1125  if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
1126    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
1127      << pivotRow_
1128      << CoinMessageEol;
1129  }
1130  // check accuracy of weights
1131  abcDualRowPivot_->checkAccuracy();
1132  // get sign for finding row of tableau
1133  // create as packed
1134#ifdef MOVE_REPLACE_PART1A
1135  if (!abcFactorization_->usingFT()) {
1136#endif
1137    usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
1138    abcFactorization_->updateColumnTranspose(usefulArray_[arrayForBtran_]);
1139#ifdef MOVE_REPLACE_PART1A
1140  } else {
1141    cilk_spawn abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
1142    usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
1143    abcFactorization_->updateColumnTransposeCpu(usefulArray_[arrayForBtran_],1);
1144    cilk_sync;
1145  }
1146#endif
1147  sequenceIn_ = -1;
1148}
1149void
1150AbcSimplexDual::getTableauColumnPart1Cilk()
1151{
1152#if ABC_DEBUG
1153  {
1154    const double * work=usefulArray_[arrayForTableauRow_].denseVector();
1155    int number=usefulArray_[arrayForTableauRow_].getNumElements();
1156    const int * which=usefulArray_[arrayForTableauRow_].getIndices();
1157    for (int i=0;i<number;i++) {
1158      if (which[i]==sequenceIn_) {
1159        assert (alpha_==work[i]);
1160        break;
1161      }
1162    }
1163  }
1164#endif
1165  //int returnCode=0;
1166  // update the incoming column
1167  unpack(usefulArray_[arrayForFtran_]);
1168  // Array[0] may be needed until updateWeights2
1169  // and update dual weights (can do in parallel - with extra array)
1170  abcFactorization_->updateColumnFTPart1(usefulArray_[arrayForFtran_]);
1171}
1172int 
1173AbcSimplexDual::getTableauColumnFlipAndStartReplaceCilk()
1174{
1175  // move checking stuff down into called functions
1176  // threads 2 and 3 are available
1177  int numberFlipped;
1178  //cilk
1179  getTableauColumnPart1Cilk();
1180#ifndef MOVE_REPLACE_PART1A
1181  cilk_spawn getTableauColumnPart2();
1182  cilk_spawn checkReplacePart1();
1183  numberFlipped=flipBounds();
1184  cilk_sync;
1185#else
1186  if (abcFactorization_->usingFT()) {
1187    cilk_spawn getTableauColumnPart2();
1188    ftAlpha_ = cilk_spawn abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
1189    numberFlipped=flipBounds();
1190    cilk_sync;
1191  } else {
1192    cilk_spawn getTableauColumnPart2();
1193    numberFlipped=flipBounds();
1194    cilk_sync;
1195  }
1196#endif
1197  //usefulArray_[arrayForTableauRow_].compact();
1198  // returns 3 if skip this iteration and re-factorize
1199  /*
1200    return code
1201    0 - OK
1202    2 - flag something and skip
1203    3 - break and re-factorize
1204    4 - break and say infeasible
1205 */
1206  int returnCode=0;
1207  // amount primal will move
1208  movement_ = -dualOut_ * directionOut_ / alpha_;
1209  // see if update stable
1210#if ABC_NORMAL_DEBUG>3
1211  if ((handler_->logLevel() & 32)) {
1212    double ft=ftAlpha_*abcFactorization_->pivotRegion()[pivotRow_];
1213    double diff1=fabs(alpha_-btranAlpha_);
1214    double diff2=fabs(fabs(alpha_)-fabs(ft));
1215    double diff3=fabs(fabs(ft)-fabs(btranAlpha_));
1216    double largest=CoinMax(CoinMax(diff1,diff2),diff3);
1217    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n", 
1218           btranAlpha_, alpha_,ft,largest);
1219    if (largest>0.001*fabs(alpha_)) {
1220      printf("bad\n");
1221    }
1222  }
1223#endif
1224  int numberPivots=abcFactorization_->pivots();
1225  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
1226  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
1227  // if can't trust much and long way from optimal then relax
1228  if (largestPrimalError_ > 10.0)
1229    checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
1230  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1231      fabs(btranAlpha_ - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
1232    handler_->message(CLP_DUAL_CHECK, messages_)
1233      << btranAlpha_
1234      << alpha_
1235      << CoinMessageEol;
1236    // see with more relaxed criterion
1237    double test;
1238    if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
1239      test = 1.0e-1 * fabs(alpha_);
1240    else
1241      test = CoinMin(10.0*checkValue,1.0e-4 * (1.0 + fabs(alpha_)));
1242    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1243                     fabs(btranAlpha_ - alpha_) > test);
1244    double derror = CoinMin(fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),1.0)*0.9999e7;
1245    int error=0;
1246    while (derror>1.0) {
1247      error++;
1248      derror *= 0.1;
1249    }
1250    int frequency[8]={99999999,100,10,2,1,1,1,1};
1251    int newFactorFrequency=CoinMin(forceFactorization_,frequency[error]);
1252#if ABC_NORMAL_DEBUG>0
1253    if (newFactorFrequency<forceFactorization_)
1254      printf("Error of %g after %d pivots old forceFact %d now %d\n",
1255             fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),numberPivots,
1256             forceFactorization_,newFactorFrequency);
1257#endif
1258    if (!numberPivots&&fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_))>0.5e-6
1259        &&abcFactorization_->pivotTolerance()<0.5)
1260      abcFactorization_->saferTolerances (1.0, -1.03);
1261    forceFactorization_=newFactorFrequency;
1262
1263   
1264#if ABC_NORMAL_DEBUG>0
1265    if (fabs(btranAlpha_ + alpha_) < 1.0e-5*(1.0 + fabs(alpha_))) {
1266      printf("diff (%g,%g) %g check %g\n",btranAlpha_,alpha_,fabs(btranAlpha_-alpha_),checkValue*(1.0+fabs(alpha_)));
1267    }
1268#endif
1269    if (numberPivots) {
1270      if (needFlag&&numberPivots<10) {
1271        // need to reject something
1272        assert (sequenceOut_>=0);
1273        char x = isColumn(sequenceOut_) ? 'C' : 'R';
1274        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1275          << x << sequenceWithin(sequenceOut_)
1276          << CoinMessageEol;
1277#if ABC_NORMAL_DEBUG>0
1278        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
1279               btranAlpha_, alpha_,numberPivots);
1280#endif
1281        // Make safer?
1282        abcFactorization_->saferTolerances (1.0, -1.03);
1283        setFlagged(sequenceOut_);
1284        // so can't happen immediately again
1285        sequenceOut_=-1;
1286        lastBadIteration_ = numberIterations_; // say be more cautious
1287      }
1288      // Make safer?
1289      //if (numberPivots<5)
1290      //abcFactorization_->saferTolerances (-0.99, -1.01);
1291      problemStatus_ = -2; // factorize now
1292      returnCode = -2;
1293      stateOfIteration_=2;
1294    } else {
1295      if (needFlag) {
1296        assert (sequenceOut_>=0);
1297        // need to reject something
1298        char x = isColumn(sequenceOut_) ? 'C' : 'R';
1299        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1300          << x << sequenceWithin(sequenceOut_)
1301          << CoinMessageEol;
1302#if ABC_NORMAL_DEBUG>3
1303        printf("flag a %g %g\n", btranAlpha_, alpha_);
1304#endif
1305        setFlagged(sequenceOut_);
1306        // so can't happen immediately again
1307        sequenceOut_=-1;
1308        //abcProgress_.clearBadTimes();
1309        lastBadIteration_ = numberIterations_; // say be more cautious
1310        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
1311          //printf("I think should declare infeasible\n");
1312          problemStatus_ = 1;
1313          returnCode = 1;
1314          stateOfIteration_=2;
1315        } else {
1316          stateOfIteration_=1;
1317        }
1318        // Make safer?
1319        if (abcFactorization_->pivotTolerance()<0.999&&stateOfIteration_==1) {
1320          // change tolerance and re-invert
1321          abcFactorization_->saferTolerances (1.0, -1.03);
1322          problemStatus_ = -6; // factorize now
1323          returnCode = -2;
1324          stateOfIteration_=2;
1325        }
1326      }
1327    }
1328  }
1329  if (!stateOfIteration_) {
1330    // check update
1331    //int updateStatus =
1332/*
1333  returns
1334  0 - OK
1335  1 - take iteration then re-factorize
1336  2 - flag something and skip
1337  3 - break and re-factorize
1338  5 - take iteration then re-factorize because of memory
1339 */
1340    int status=checkReplace();
1341    if (status&&!returnCode)
1342      returnCode=status;
1343  }
1344  //clearArrays(arrayForFlipRhs_);
1345  if (stateOfIteration_&&numberFlipped) {
1346    //usefulArray_[arrayForTableauRow_].compact();
1347    // move variables back
1348    flipBack(numberFlipped);
1349  }
1350  // could choose average of all three
1351  return returnCode;
1352}
1353int 
1354AbcSimplexDual::whileIteratingParallel(int numberIterations)
1355{
1356  int returnCode=-1;
1357#ifdef EARLY_FACTORIZE
1358  int savePivot=-1;
1359  CoinIndexedVector & early = usefulArray_[ABC_NUMBER_USEFUL-1];
1360  int * pivotIndices = early.getIndices();
1361  double * pivotValues = early.denseVector();
1362  int pivotCountPosition=early.capacity()-1;
1363  int numberSaved=0;
1364  if (numberIterations==123456789)
1365    savePivot=pivotCountPosition;;
1366#endif
1367  numberIterations += numberIterations_;
1368  do {
1369#if ABC_DEBUG
1370    checkArrays();
1371#endif
1372    /*
1373      -1 - just move dual in values pass
1374      0 - take iteration
1375      1 - don't take but continue
1376      2 - don't take and break
1377    */
1378    stateOfIteration_=0;
1379    returnCode=-1;
1380    // put row of tableau in usefulArray[arrayForTableauRow_]
1381    /*
1382      Could
1383      a) copy [arrayForBtran_] and start off updateWeights earlier
1384      b) break transposeTimes into two and do after slack part
1385      c) also think if cleaner to go all way with update but just don't do final part
1386    */
1387    //upperTheta=COIN_DBL_MAX;
1388    double saveAcceptable=acceptablePivot_;
1389    if (largestPrimalError_>1.0e-5||largestDualError_>1.0e-5) {
1390    //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) {
1391      if (!abcFactorization_->pivots())
1392        acceptablePivot_*=1.0e2;
1393      else if (abcFactorization_->pivots()<5)
1394        acceptablePivot_*=1.0e1;
1395    }
1396#ifdef MOVE_UPDATE_WEIGHTS
1397    // copy btran across
1398    usefulArray_[5].copy(usefulArray_[arrayForBtran_]);
1399    cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[5]);;
1400#endif
1401    dualColumn1();
1402    acceptablePivot_=saveAcceptable;
1403    if ((stateOfProblem_&VALUES_PASS)!=0) {
1404      // see if can just move dual
1405      if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
1406        stateOfIteration_=-1;
1407      }
1408    }
1409    if (!stateOfIteration_) {
1410#ifndef MOVE_UPDATE_WEIGHTS
1411      cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[arrayForBtran_]);;
1412#endif
1413      // get sequenceIn_
1414      dualPivotColumn();
1415      //stateOfIteration_=0;
1416      if (sequenceIn_ >= 0) {
1417        // normal iteration
1418        // update the incoming column
1419        //arrayForReplaceColumn_=getAvailableArray();
1420        getTableauColumnFlipAndStartReplaceCilk();
1421        //usleep(1000);
1422      } else if (stateOfIteration_!=-1) {
1423        stateOfIteration_=2;
1424      }
1425    }
1426    cilk_sync;
1427    // Check event
1428    {
1429      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1430      if (status >= 0) {
1431        problemStatus_ = 5;
1432        secondaryStatus_ = ClpEventHandler::endOfIteration;
1433        returnCode = 4;
1434        stateOfIteration_=2;
1435      }
1436    }
1437    if (!stateOfIteration_) {
1438      //        assert (stateOfIteration_!=1);
1439      int whatNext=0;
1440      whatNext = housekeeping();
1441#ifdef EARLY_FACTORIZE
1442      if (savePivot>=0) {
1443        // save pivot
1444        pivotIndices[--savePivot]=sequenceOut_;
1445        pivotValues[savePivot]=alpha_;
1446        pivotIndices[--savePivot]=sequenceIn_;
1447        pivotValues[savePivot]=btranAlpha_;
1448        numberSaved++;
1449      }
1450#endif
1451      if (whatNext == 1) {
1452        problemStatus_ = -2; // refactorize
1453        usefulArray_[arrayForTableauRow_].compact();
1454      } else if (whatNext == 2) {
1455        // maximum iterations or equivalent
1456        problemStatus_ = 3;
1457        returnCode = 3;
1458        stateOfIteration_=2;
1459      }
1460#ifdef EARLY_FACTORIZE
1461      if (savePivot>=0) {
1462        // tell worker can update this
1463        pivotIndices[pivotCountPosition]=numberSaved;
1464      }
1465#endif
1466    } else {
1467      usefulArray_[arrayForTableauRow_].compact();
1468       if (stateOfIteration_<0) {
1469        // move dual in dual values pass
1470        theta_=abcDj_[sequenceOut_];
1471        updateDualsInDual();
1472        abcDj_[sequenceOut_]=0.0;
1473        sequenceOut_=-1;
1474      }
1475      // clear all arrays
1476      clearArrays(-1);
1477    }
1478    // at this stage sequenceIn_ may be <0
1479    if (sequenceIn_<0&&sequenceOut_>=0) {
1480      usefulArray_[arrayForTableauRow_].compact();
1481      // no incoming column is valid
1482      returnCode=noPivotColumn();
1483    }
1484    if (stateOfIteration_==2) {
1485      usefulArray_[arrayForTableauRow_].compact();
1486      sequenceIn_=-1;
1487#ifdef ABC_LONG_FACTORIZATION
1488      abcFactorization_->clearHiddenArrays();
1489#endif
1490      break;
1491    }
1492    if (problemStatus_!=-1) {
1493#ifndef MOVE_UPDATE_WEIGHTS
1494      abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
1495#else
1496      abcDualRowPivot_->updateWeights2(usefulArray_[5],usefulArray_[arrayForFtran_]);
1497#endif
1498#ifdef ABC_LONG_FACTORIZATION
1499      abcFactorization_->clearHiddenArrays();
1500#endif
1501      swapPrimalStuff();
1502      break;
1503    }
1504    if (stateOfIteration_==1) {
1505      // clear all arrays
1506      clearArrays(-2);
1507#ifdef ABC_LONG_FACTORIZATION
1508      abcFactorization_->clearHiddenArrays();
1509#endif
1510    }
1511    if (stateOfIteration_==0) {
1512      // can do these in parallel
1513      // No idea why I need this - but otherwise runs not repeatable (try again??)
1514      //usefulArray_[3].compact();
1515      cilk_spawn updateDualsInDual();
1516      int lastSequenceOut;
1517      int lastDirectionOut;
1518      if (firstFree_<0) {
1519        // can do in parallel
1520        cilk_spawn replaceColumnPart3();
1521        updatePrimalSolution();
1522        swapPrimalStuff();
1523        // dualRow will go to virtual row pivot choice algorithm
1524        // make sure values pass off if it should be
1525        // use Btran array and clear inside dualPivotRow (if used)
1526        lastSequenceOut=sequenceOut_;
1527        lastDirectionOut=directionOut_;
1528        dualPivotRow();
1529        cilk_sync;
1530      } else {
1531        // be more careful as dualPivotRow may do update
1532        cilk_spawn replaceColumnPart3();
1533        updatePrimalSolution();
1534        swapPrimalStuff();
1535        // dualRow will go to virtual row pivot choice algorithm
1536        // make sure values pass off if it should be
1537        // use Btran array and clear inside dualPivotRow (if used)
1538        lastSequenceOut=sequenceOut_;
1539        lastDirectionOut=directionOut_;
1540        cilk_sync;
1541        dualPivotRow();
1542      }
1543      lastPivotRow_=pivotRow_;
1544      if (pivotRow_ >= 0) {
1545        // we found a pivot row
1546        // create as packed
1547        // dual->checkReplacePart1a();
1548        createDualPricingVectorCilk();
1549        swapDualStuff(lastSequenceOut,lastDirectionOut);
1550      }
1551      cilk_sync;
1552    } else {
1553      // after moving dual in values pass
1554      dualPivotRow();
1555      lastPivotRow_=pivotRow_;
1556      if (pivotRow_ >= 0) {
1557        // we found a pivot row
1558        // create as packed
1559        createDualPricingVectorCilk();
1560      }
1561    }
1562    if (pivotRow_<0) {
1563      // no pivot row
1564      clearArrays(-1);
1565      returnCode=noPivotRow();
1566      break;
1567    }
1568    if (numberIterations == numberIterations_&&problemStatus_==-1) {
1569      returnCode=-99;
1570      break;
1571    }
1572  } while (problemStatus_ == -1);
1573  return returnCode;
1574}
1575#if 0
1576void
1577AbcSimplexDual::whileIterating2()
1578{
1579  // get sequenceIn_
1580  dualPivotColumn();
1581  //stateOfIteration_=0;
1582  if (sequenceIn_ >= 0) {
1583    // normal iteration
1584    // update the incoming column
1585    //arrayForReplaceColumn_=getAvailableArray();
1586    getTableauColumnFlipAndStartReplaceCilk();
1587    //usleep(1000);
1588  } else if (stateOfIteration_!=-1) {
1589    stateOfIteration_=2;
1590  }
1591}
1592#endif
1593#endif
1594void 
1595AbcSimplexDual::updatePrimalSolution()
1596{
1597  // finish doing weights
1598#ifndef MOVE_UPDATE_WEIGHTS
1599  abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
1600                                         movement_);
1601#else
1602  abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[5],usefulArray_[arrayForFtran_],
1603                                         movement_);
1604#endif
1605}
1606int xxInfo[6][8];
1607double parallelDual4(AbcSimplexDual * dual)
1608{
1609  int maximumRows=dual->maximumAbcNumberRows();
1610  //int numberTotal=dual->numberTotal();
1611  CoinIndexedVector & update = *dual->usefulArray(dual->arrayForBtran());
1612  CoinPartitionedVector & tableauRow= *dual->usefulArray(dual->arrayForTableauRow());
1613  CoinPartitionedVector & candidateList= *dual->usefulArray(dual->arrayForDualColumn());
1614  AbcMatrix * matrix = dual->abcMatrix();
1615  // do slacks first (move before sync)
1616  double upperTheta=dual->upperTheta();
1617  //assert (upperTheta>-100*dual->dualTolerance()||dual->sequenceIn()>=0||dual->lastFirstFree()>=0);
1618#if ABC_PARALLEL
1619#define NUMBER_BLOCKS NUMBER_ROW_BLOCKS
1620  int numberBlocks=CoinMin(NUMBER_BLOCKS,dual->numberCpus());
1621#else
1622#define NUMBER_BLOCKS 1
1623  int numberBlocks=NUMBER_BLOCKS;
1624#endif
1625#if ABC_PARALLEL
1626  if (dual->parallelMode()==0)
1627    numberBlocks=1;
1628#endif
1629  // see if worth using row copy
1630  bool gotRowCopy = matrix->gotRowCopy();
1631  int number=update.getNumElements();
1632  assert (number);
1633  bool useRowCopy = (gotRowCopy&&(number<<2)<maximumRows);
1634  assert (numberBlocks==matrix->numberRowBlocks());
1635#if ABC_PARALLEL==1
1636  // redo so can pass info in with void *
1637  CoinAbcThreadInfo * infoP = dual->threadInfoPointer();
1638  int cpuMask=((1<<dual->parallelMode())-1);
1639  cpuMask += cpuMask<<5;
1640  dual->setStopStart(cpuMask);
1641#endif
1642  CoinAbcThreadInfo info[NUMBER_BLOCKS];
1643  const int * starts;
1644  if (useRowCopy)
1645    starts=matrix->blockStart();
1646  else
1647    starts=matrix->startColumnBlock();
1648  tableauRow.setPartitions(numberBlocks,starts);
1649  candidateList.setPartitions(numberBlocks,starts);
1650  //printf("free sequence in %d\n",dual->freeSequenceIn());
1651  if (useRowCopy) {
1652    // using row copy
1653#if ABC_PARALLEL
1654    if (numberBlocks>1) {
1655#if ABC_PARALLEL==2
1656    for (int i=0;i<numberBlocks;i++) {
1657      info[i].stuff[1]=i;
1658      info[i].stuff[2]=-1;
1659      info[i].result=upperTheta;
1660      info[i].result= cilk_spawn
1661        matrix->dualColumn1Row(info[i].stuff[1],COIN_DBL_MAX,info[i].stuff[2],
1662                               update,tableauRow,candidateList);
1663    }
1664    cilk_sync;
1665#else
1666      // parallel 1
1667    for (int i=0;i<numberBlocks;i++) {
1668      info[i].status=5;
1669      info[i].stuff[1]=i;
1670      info[i].result=upperTheta;
1671      if (i<numberBlocks-1) {
1672        infoP[i]=info[i];
1673        if (i==numberBlocks-2) 
1674          dual->startParallelStuff(5);
1675      }
1676    }
1677    info[numberBlocks-1].stuff[1]=numberBlocks-1;
1678    info[numberBlocks-1].stuff[2]=-1;
1679    //double freeAlpha;
1680    info[numberBlocks-1].result =
1681      matrix->dualColumn1Row(info[numberBlocks-1].stuff[1],
1682                             upperTheta,info[numberBlocks-1].stuff[2],
1683                             update,tableauRow,candidateList);
1684    dual->stopParallelStuff(5);
1685    for (int i=0;i<numberBlocks-1;i++) 
1686      info[i]=infoP[i];
1687#endif
1688    } else {
1689#endif
1690      info[0].status=5;
1691      info[0].stuff[1]=0;
1692      info[0].result=upperTheta;
1693      info[0].stuff[1]=0;
1694      info[0].stuff[2]=-1;
1695      // worth using row copy
1696      //assert (number>2);
1697      info[0].result =
1698      matrix->dualColumn1Row(info[0].stuff[1],
1699                             upperTheta,info[0].stuff[2],
1700                             update,tableauRow,candidateList);
1701#if ABC_PARALLEL
1702    }
1703#endif
1704  } else { // end of useRowCopy
1705#if ABC_PARALLEL
1706    if (numberBlocks>1) {
1707#if ABC_PARALLEL==2
1708      // do by column
1709      for (int i=0;i<numberBlocks;i++) {
1710        info[i].stuff[1]=i;
1711        info[i].result=upperTheta;
1712        cilk_spawn
1713          matrix->dualColumn1Part(info[i].stuff[1],info[i].stuff[2],
1714                                             info[i].result,
1715                                             update,tableauRow,candidateList);
1716      }
1717      cilk_sync;
1718#else
1719      // parallel 1
1720      // do by column
1721      for (int i=0;i<numberBlocks;i++) {
1722        info[i].status=4;
1723        info[i].stuff[1]=i;
1724        info[i].result=upperTheta;
1725        if (i<numberBlocks-1) {
1726          infoP[i]=info[i];
1727          if (i==numberBlocks-2) 
1728            dual->startParallelStuff(4);
1729        }
1730      }
1731      matrix->dualColumn1Part(info[numberBlocks-1].stuff[1],info[numberBlocks-1].stuff[2],
1732                                                            info[numberBlocks-1].result,
1733                                                            update,tableauRow,candidateList);
1734      dual->stopParallelStuff(4);
1735      for (int i=0;i<numberBlocks-1;i++) 
1736        info[i]=infoP[i];
1737#endif
1738    } else {
1739#endif
1740      info[0].status=4;
1741      info[0].stuff[1]=0;
1742      info[0].result=upperTheta;
1743      info[0].stuff[1]=0;
1744      matrix->dualColumn1Part(info[0].stuff[1],info[0].stuff[2],
1745                                               info[0].result,
1746                                               update,tableauRow,candidateList);
1747#if ABC_PARALLEL
1748    }
1749#endif
1750  }
1751  int sequenceIn[NUMBER_BLOCKS];
1752  bool anyFree=false;
1753#if 0
1754  if (useRowCopy) {
1755    printf("Result at iteration %d slack seqIn %d upperTheta %g\n",
1756           dual->numberIterations(),dual->freeSequenceIn(),upperTheta);
1757    double * arrayT = tableauRow.denseVector();
1758    int * indexT = tableauRow.getIndices();
1759    double * arrayC = candidateList.denseVector();
1760    int * indexC = candidateList.getIndices();
1761    for (int i=0;i<numberBlocks;i++) {
1762      printf("Block %d seqIn %d upperTheta %g first %d last %d firstIn %d nnz %d numrem %d\n",
1763             i,info[i].stuff[2],info[i].result,
1764             xxInfo[0][i],xxInfo[1][i],xxInfo[2][i],xxInfo[3][i],xxInfo[4][i]);
1765      if (xxInfo[3][i]<-35) {
1766        assert (xxInfo[3][i]==tableauRow.getNumElements(i));
1767        assert (xxInfo[4][i]==candidateList.getNumElements(i));
1768        for (int k=0;k<xxInfo[3][i];k++)
1769          printf("T %d %d %g\n",k,indexT[k+xxInfo[2][i]],arrayT[k+xxInfo[2][i]]);
1770        for (int k=0;k<xxInfo[4][i];k++)
1771          printf("C %d %d %g\n",k,indexC[k+xxInfo[2][i]],arrayC[k+xxInfo[2][i]]);
1772      }
1773    }
1774  }
1775#endif
1776  for (int i=0;i<numberBlocks;i++) {
1777    sequenceIn[i]=info[i].stuff[2];
1778    if (sequenceIn[i]>=0)
1779      anyFree=true;
1780    upperTheta=CoinMin(info[i].result,upperTheta);
1781    //assert (info[i].result>-100*dual->dualTolerance()||sequenceIn[i]>=0||dual->lastFirstFree()>=0);
1782  }
1783  if (anyFree) {
1784    int *  COIN_RESTRICT index = tableauRow.getIndices();
1785    double *  COIN_RESTRICT array = tableauRow.denseVector();
1786    // search for free coming in
1787    double bestFree=0.0;
1788    int bestSequence=dual->sequenceIn();
1789    if (bestSequence>=0)
1790      bestFree=dual->alpha();
1791    for (int i=0;i<numberBlocks;i++) {
1792      int iLook=sequenceIn[i];
1793      if (iLook>=0) {
1794        // free variable - search
1795        int start=starts[i];
1796        int end=start+tableauRow.getNumElements(i);
1797        for (int k=start;k<end;k++) {
1798          if (iLook==index[k]) {
1799            if (fabs(bestFree)<fabs(array[k])) {
1800              bestFree=array[k];
1801              bestSequence=iLook;
1802            }
1803            break;
1804          }
1805        }
1806      }
1807    }
1808    if (bestSequence>=0) {
1809      double oldValue = dual->djRegion()[bestSequence];
1810      dual->setSequenceIn(bestSequence);
1811      dual->setAlpha(bestFree);
1812      dual->setTheta(oldValue / bestFree);
1813    }
1814  }
1815  //tableauRow.compact();
1816  //candidateList.compact();
1817#if 0//ndef NDEBUG
1818  tableauRow.setPackedMode(true);
1819  tableauRow.sortPacked();
1820  candidateList.setPackedMode(true);
1821  candidateList.sortPacked();
1822#endif
1823  return upperTheta;
1824}
1825#if ABC_PARALLEL==2
1826static void
1827parallelDual5a(AbcSimplexFactorization * factorization,
1828              CoinIndexedVector * whichVector,
1829              int numberCpu,
1830              int whichCpu,
1831              double * weights)
1832{
1833  double *  COIN_RESTRICT array = whichVector->denseVector();
1834  int *  COIN_RESTRICT which = whichVector->getIndices();
1835  int numberRows=factorization->numberRows();
1836  for (int i = whichCpu; i < numberRows; i+=numberCpu) {
1837    double value = 0.0;
1838    array[i] = 1.0;
1839    which[0] = i;
1840    whichVector->setNumElements(1);
1841    factorization->updateColumnTransposeCpu(*whichVector,whichCpu);
1842    int number = whichVector->getNumElements();
1843    for (int j = 0; j < number; j++) {
1844      int k= which[j];
1845      value += array[k] * array[k];
1846      array[k] = 0.0;
1847    }
1848    weights[i] = value;
1849  }
1850  whichVector->setNumElements(0);
1851}
1852#endif
1853#if ABC_PARALLEL==2
1854void
1855parallelDual5(AbcSimplexFactorization * factorization,
1856              CoinIndexedVector ** whichVector,
1857              int numberCpu,
1858              int whichCpu,
1859              double * weights)
1860{
1861  if (whichCpu) {
1862    cilk_spawn parallelDual5(factorization,whichVector,numberCpu,whichCpu-1,weights);
1863    parallelDual5a(factorization,whichVector[whichCpu],numberCpu,whichCpu,weights);
1864    cilk_sync;
1865  } else {
1866    parallelDual5a(factorization,whichVector[whichCpu],numberCpu,whichCpu,weights);
1867  }
1868}
1869#endif
1870// cilk seems a bit fragile
1871#define CILK_FRAGILE 1
1872#if CILK_FRAGILE>1
1873#undef cilk_spawn
1874#undef cilk_sync
1875#define cilk_spawn
1876#define cilk_sync
1877#define ONWARD 0
1878#elif CILK_FRAGILE==1
1879#define ONWARD 0
1880#else
1881#define ONWARD 1
1882#endif
1883// Code below is just a translation of LAPACK
1884#define BLOCKING1 8 // factorization strip
1885#define BLOCKING2 8 // dgemm recursive
1886#define BLOCKING3 16 // dgemm parallel
1887/* type
1888   0 Left Lower NoTranspose Unit
1889   1 Left Upper NoTranspose NonUnit
1890   2 Left Lower Transpose Unit
1891   3 Left Upper Transpose NonUnit
1892*/
1893static void CoinAbcDtrsmFactor(int m, int n, double * COIN_RESTRICT a,int lda)
1894{
1895  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
1896  assert (m==BLOCKING8);
1897  // 0 Left Lower NoTranspose Unit
1898  /* entry for column j and row i (when multiple of BLOCKING8)
1899     is at aBlocked+j*m+i*BLOCKING8
1900  */
1901  double * COIN_RESTRICT aBase2 = a;
1902  double * COIN_RESTRICT bBase2 = aBase2+lda*BLOCKING8;
1903  for (int jj=0;jj<n;jj+=BLOCKING8) {
1904    double * COIN_RESTRICT bBase = bBase2;
1905    for (int j=jj;j<jj+BLOCKING8;j++) {
1906#if 0
1907      double * COIN_RESTRICT aBase = aBase2;
1908      for (int k=0;k<BLOCKING8;k++) {
1909        double bValue = bBase[k];
1910        if (bValue) {
1911          for (int i=k+1;i<BLOCKING8;i++) {
1912            bBase[i]-=bValue*aBase[i];
1913          }
1914        }
1915        aBase+=BLOCKING8;
1916      }
1917#else
1918      // unroll - stay in registers - don't test for zeros
1919      assert (BLOCKING8==8);
1920      double bValue0 = bBase[0];
1921      double bValue1 = bBase[1];
1922      double bValue2 = bBase[2];
1923      double bValue3 = bBase[3];
1924      double bValue4 = bBase[4];
1925      double bValue5 = bBase[5];
1926      double bValue6 = bBase[6];
1927      double bValue7 = bBase[7];
1928      bValue1-=bValue0*a[1+0*BLOCKING8];
1929      bBase[1] = bValue1;
1930      bValue2-=bValue0*a[2+0*BLOCKING8];
1931      bValue3-=bValue0*a[3+0*BLOCKING8];
1932      bValue4-=bValue0*a[4+0*BLOCKING8];
1933      bValue5-=bValue0*a[5+0*BLOCKING8];
1934      bValue6-=bValue0*a[6+0*BLOCKING8];
1935      bValue7-=bValue0*a[7+0*BLOCKING8];
1936      bValue2-=bValue1*a[2+1*BLOCKING8];
1937      bBase[2] = bValue2;
1938      bValue3-=bValue1*a[3+1*BLOCKING8];
1939      bValue4-=bValue1*a[4+1*BLOCKING8];
1940      bValue5-=bValue1*a[5+1*BLOCKING8];
1941      bValue6-=bValue1*a[6+1*BLOCKING8];
1942      bValue7-=bValue1*a[7+1*BLOCKING8];
1943      bValue3-=bValue2*a[3+2*BLOCKING8];
1944      bBase[3] = bValue3;
1945      bValue4-=bValue2*a[4+2*BLOCKING8];
1946      bValue5-=bValue2*a[5+2*BLOCKING8];
1947      bValue6-=bValue2*a[6+2*BLOCKING8];
1948      bValue7-=bValue2*a[7+2*BLOCKING8];
1949      bValue4-=bValue3*a[4+3*BLOCKING8];
1950      bBase[4] = bValue4;
1951      bValue5-=bValue3*a[5+3*BLOCKING8];
1952      bValue6-=bValue3*a[6+3*BLOCKING8];
1953      bValue7-=bValue3*a[7+3*BLOCKING8];
1954      bValue5-=bValue4*a[5+4*BLOCKING8];
1955      bBase[5] = bValue5;
1956      bValue6-=bValue4*a[6+4*BLOCKING8];
1957      bValue7-=bValue4*a[7+4*BLOCKING8];
1958      bValue6-=bValue5*a[6+5*BLOCKING8];
1959      bBase[6] = bValue6;
1960      bValue7-=bValue5*a[7+5*BLOCKING8];
1961      bValue7-=bValue6*a[7+6*BLOCKING8];
1962      bBase[7] = bValue7;
1963#endif
1964      bBase += BLOCKING8;
1965    }
1966    bBase2 +=lda*BLOCKING8; 
1967  }
1968}
1969#define UNROLL_DTRSM 16
1970#define CILK_DTRSM 32
1971static void dtrsm0(int kkk, int first, int last,
1972                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
1973{
1974  int mm=CoinMin(kkk+UNROLL_DTRSM*BLOCKING8,m);
1975  assert ((last-first)%BLOCKING8==0);
1976  if (last-first>CILK_DTRSM) {
1977    int mid=((first+last)>>4)<<3;
1978    cilk_spawn dtrsm0(kkk,first,mid,m,a,b);
1979    dtrsm0(kkk,mid,last,m,a,b);
1980    cilk_sync;
1981  } else {
1982    const double * COIN_RESTRICT aBaseA = a+UNROLL_DTRSM*BLOCKING8X8+kkk*BLOCKING8;
1983    aBaseA += (first-mm)*BLOCKING8-BLOCKING8X8;
1984    aBaseA += m*kkk;
1985    for (int ii=first;ii<last;ii+=BLOCKING8) {
1986      aBaseA += BLOCKING8X8;
1987      const double * COIN_RESTRICT aBaseB=aBaseA;
1988      double * COIN_RESTRICT bBaseI = b+ii;
1989      for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
1990        double * COIN_RESTRICT bBase = b+kk;
1991        const double * COIN_RESTRICT aBase2 = aBaseB;//a+UNROLL_DTRSM*BLOCKING8X8+m*kk+kkk*BLOCKING8;
1992        //aBase2 += (ii-mm)*BLOCKING8;
1993        //assert (aBase2==aBaseB);
1994        aBaseB+=m*BLOCKING8;
1995#if AVX2 !=2
1996#define ALTERNATE_INNER
1997#ifndef ALTERNATE_INNER
1998        for (int k=0;k<BLOCKING8;k++) {
1999          //coin_prefetch_const(aBase2+BLOCKING8);             
2000          for (int i=0;i<BLOCKING8;i++) {
2001            bBaseI[i] -= bBase[k]*aBase2[i];
2002          }
2003          aBase2 += BLOCKING8;
2004        }
2005#else
2006        double b0=bBaseI[0];
2007        double b1=bBaseI[1];
2008        double b2=bBaseI[2];
2009        double b3=bBaseI[3];
2010        double b4=bBaseI[4];
2011        double b5=bBaseI[5];
2012        double b6=bBaseI[6];
2013        double b7=bBaseI[7];
2014        for (int k=0;k<BLOCKING8;k++) {
2015          //coin_prefetch_const(aBase2+BLOCKING8);             
2016          double bValue=bBase[k];
2017          b0 -= bValue * aBase2[0];
2018          b1 -= bValue * aBase2[1];
2019          b2 -= bValue * aBase2[2];
2020          b3 -= bValue * aBase2[3];
2021          b4 -= bValue * aBase2[4];
2022          b5 -= bValue * aBase2[5];
2023          b6 -= bValue * aBase2[6];
2024          b7 -= bValue * aBase2[7];
2025          aBase2 += BLOCKING8;
2026        }
2027        bBaseI[0]=b0;
2028        bBaseI[1]=b1;
2029        bBaseI[2]=b2;
2030        bBaseI[3]=b3;
2031        bBaseI[4]=b4;
2032        bBaseI[5]=b5;
2033        bBaseI[6]=b6;
2034        bBaseI[7]=b7;
2035#endif
2036#else
2037        __m256d b0=_mm256_load_pd(bBaseI);
2038        __m256d b1=_mm256_load_pd(bBaseI+4);
2039        for (int k=0;k<BLOCKING8;k++) {
2040          //coin_prefetch_const(aBase2+BLOCKING8);             
2041          __m256d bb = _mm256_broadcast_sd(bBase+k);
2042          __m256d a0 = _mm256_load_pd(aBase2);
2043          b0 -= a0*bb;
2044          __m256d a1 = _mm256_load_pd(aBase2+4);
2045          b1 -= a1*bb;
2046          aBase2 += BLOCKING8;
2047        }
2048        //_mm256_store_pd ((bBaseI), (b0));
2049        *reinterpret_cast<__m256d *>(bBaseI)=b0;
2050        //_mm256_store_pd (bBaseI+4, b1);
2051        *reinterpret_cast<__m256d *>(bBaseI+4)=b1;
2052#endif
2053      }
2054    }
2055  }
2056}
2057void CoinAbcDtrsm0(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2058{
2059  assert ((m&(BLOCKING8-1))==0);
2060  // 0 Left Lower NoTranspose Unit
2061  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM*BLOCKING8) {
2062    int mm=CoinMin(m,kkk+UNROLL_DTRSM*BLOCKING8);
2063    for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
2064      const double * COIN_RESTRICT aBase2 = a+kk*(m+BLOCKING8);
2065      double * COIN_RESTRICT bBase = b+kk;
2066      for (int k=0;k<BLOCKING8;k++) {
2067        for (int i=k+1;i<BLOCKING8;i++) {
2068          bBase[i] -= bBase[k]*aBase2[i];
2069        }
2070        aBase2 += BLOCKING8;
2071      }
2072      for (int ii=kk+BLOCKING8;ii<mm;ii+=BLOCKING8) {
2073        double * COIN_RESTRICT bBaseI = b+ii;
2074        for (int k=0;k<BLOCKING8;k++) {
2075          //coin_prefetch_const(aBase2+BLOCKING8);             
2076          for (int i=0;i<BLOCKING8;i++) {
2077            bBaseI[i] -= bBase[k]*aBase2[i];
2078          }
2079          aBase2 += BLOCKING8;
2080        }
2081      }
2082    }
2083    dtrsm0(kkk,mm,m,m,a,b);
2084  }
2085}
2086static void dtrsm1(int kkk, int first, int last,
2087                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2088{
2089  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2090  assert ((last-first)%BLOCKING8==0);
2091  if (last-first>CILK_DTRSM) {
2092    int mid=((first+last)>>4)<<3;
2093    cilk_spawn dtrsm1(kkk,first,mid,m,a,b);
2094    dtrsm1(kkk,mid,last,m,a,b);
2095    cilk_sync;
2096  } else {
2097    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2098      double * COIN_RESTRICT bBase2 = b+iii;
2099      const double * COIN_RESTRICT aBaseA=
2100        a+BLOCKING8X8+BLOCKING8*iii;
2101      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2102        double * COIN_RESTRICT bBase = b+ii;
2103        const double * COIN_RESTRICT aBase=aBaseA+m*ii;
2104#if AVX2 !=2
2105#ifndef ALTERNATE_INNER
2106        for (int i=BLOCKING8-1;i>=0;i--) {
2107          aBase -= BLOCKING8;
2108          //coin_prefetch_const(aBase-BLOCKING8);               
2109          for (int k=BLOCKING8-1;k>=0;k--) {
2110            bBase2[k] -= bBase[i]*aBase[k];
2111          }
2112        }
2113#else
2114        double b0=bBase2[0];
2115        double b1=bBase2[1];
2116        double b2=bBase2[2];
2117        double b3=bBase2[3];
2118        double b4=bBase2[4];
2119        double b5=bBase2[5];
2120        double b6=bBase2[6];
2121        double b7=bBase2[7];
2122        for (int i=BLOCKING8-1;i>=0;i--) {
2123          aBase -= BLOCKING8;
2124          //coin_prefetch_const(aBase-BLOCKING8);               
2125          double bValue=bBase[i];
2126          b0 -= bValue * aBase[0];
2127          b1 -= bValue * aBase[1];
2128          b2 -= bValue * aBase[2];
2129          b3 -= bValue * aBase[3];
2130          b4 -= bValue * aBase[4];
2131          b5 -= bValue * aBase[5];
2132          b6 -= bValue * aBase[6];
2133          b7 -= bValue * aBase[7];
2134        }
2135        bBase2[0]=b0;
2136        bBase2[1]=b1;
2137        bBase2[2]=b2;
2138        bBase2[3]=b3;
2139        bBase2[4]=b4;
2140        bBase2[5]=b5;
2141        bBase2[6]=b6;
2142        bBase2[7]=b7;
2143#endif
2144#else
2145        __m256d b0=_mm256_load_pd(bBase2);
2146        __m256d b1=_mm256_load_pd(bBase2+4);
2147        for (int i=BLOCKING8-1;i>=0;i--) {
2148          aBase -= BLOCKING8;
2149          //coin_prefetch_const(aBase5-BLOCKING8);             
2150          __m256d bb = _mm256_broadcast_sd(bBase+i);
2151          __m256d a0 = _mm256_load_pd(aBase);
2152          b0 -= a0*bb;
2153          __m256d a1 = _mm256_load_pd(aBase+4);
2154          b1 -= a1*bb;
2155        }
2156        //_mm256_store_pd (bBase2, b0);
2157        *reinterpret_cast<__m256d *>(bBase2)=b0;
2158        //_mm256_store_pd (bBase2+4, b1);
2159        *reinterpret_cast<__m256d *>(bBase2+4)=b1;
2160#endif
2161      }
2162    }
2163  }
2164}
2165void CoinAbcDtrsm1(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2166{
2167  assert ((m&(BLOCKING8-1))==0);
2168  // 1 Left Upper NoTranspose NonUnit
2169  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2170    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2171    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2172      const double * COIN_RESTRICT aBase = a+m*ii+ii*BLOCKING8;
2173      double * COIN_RESTRICT bBase = b+ii;
2174      for (int i=BLOCKING8-1;i>=0;i--) {
2175        bBase[i] *= aBase[i*(BLOCKING8+1)];
2176        for (int k=i-1;k>=0;k--) {
2177          bBase[k] -= bBase[i]*aBase[k+i*BLOCKING8];
2178        }
2179      }
2180      double * COIN_RESTRICT bBase2 = bBase;
2181      for (int iii=ii;iii>mm;iii-=BLOCKING8) {
2182        bBase2 -= BLOCKING8;
2183        for (int i=BLOCKING8-1;i>=0;i--) {
2184          aBase -= BLOCKING8;
2185          coin_prefetch_const(aBase-BLOCKING8);         
2186          for (int k=BLOCKING8-1;k>=0;k--) {
2187            bBase2[k] -= bBase[i]*aBase[k];
2188          }
2189        }
2190      }
2191    }
2192    dtrsm1(kkk, 0, mm,  m, a, b);
2193  }
2194}
2195static void dtrsm2(int kkk, int first, int last,
2196                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2197{
2198  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2199  assert ((last-first)%BLOCKING8==0);
2200  if (last-first>CILK_DTRSM) {
2201    int mid=((first+last)>>4)<<3;
2202    cilk_spawn dtrsm2(kkk,first,mid,m,a,b);
2203    dtrsm2(kkk,mid,last,m,a,b);
2204    cilk_sync;
2205  } else {
2206    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2207      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2208        const double * COIN_RESTRICT bBase = b+ii;
2209        double * COIN_RESTRICT bBase2=b+iii;
2210        const double * COIN_RESTRICT aBase=a+ii*BLOCKING8+m*iii+BLOCKING8X8;
2211        for (int j=BLOCKING8-1;j>=0;j-=4) {
2212          double bValue0=bBase2[j];
2213          double bValue1=bBase2[j-1];
2214          double bValue2=bBase2[j-2];
2215          double bValue3=bBase2[j-3];
2216          bValue0 -= aBase[-1*BLOCKING8+7]*bBase[7];
2217          bValue1 -= aBase[-2*BLOCKING8+7]*bBase[7];
2218          bValue2 -= aBase[-3*BLOCKING8+7]*bBase[7];
2219          bValue3 -= aBase[-4*BLOCKING8+7]*bBase[7];
2220          bValue0 -= aBase[-1*BLOCKING8+6]*bBase[6];
2221          bValue1 -= aBase[-2*BLOCKING8+6]*bBase[6];
2222          bValue2 -= aBase[-3*BLOCKING8+6]*bBase[6];
2223          bValue3 -= aBase[-4*BLOCKING8+6]*bBase[6];
2224          bValue0 -= aBase[-1*BLOCKING8+5]*bBase[5];
2225          bValue1 -= aBase[-2*BLOCKING8+5]*bBase[5];
2226          bValue2 -= aBase[-3*BLOCKING8+5]*bBase[5];
2227          bValue3 -= aBase[-4*BLOCKING8+5]*bBase[5];
2228          bValue0 -= aBase[-1*BLOCKING8+4]*bBase[4];
2229          bValue1 -= aBase[-2*BLOCKING8+4]*bBase[4];
2230          bValue2 -= aBase[-3*BLOCKING8+4]*bBase[4];
2231          bValue3 -= aBase[-4*BLOCKING8+4]*bBase[4];
2232          bValue0 -= aBase[-1*BLOCKING8+3]*bBase[3];
2233          bValue1 -= aBase[-2*BLOCKING8+3]*bBase[3];
2234          bValue2 -= aBase[-3*BLOCKING8+3]*bBase[3];
2235          bValue3 -= aBase[-4*BLOCKING8+3]*bBase[3];
2236          bValue0 -= aBase[-1*BLOCKING8+2]*bBase[2];
2237          bValue1 -= aBase[-2*BLOCKING8+2]*bBase[2];
2238          bValue2 -= aBase[-3*BLOCKING8+2]*bBase[2];
2239          bValue3 -= aBase[-4*BLOCKING8+2]*bBase[2];
2240          bValue0 -= aBase[-1*BLOCKING8+1]*bBase[1];
2241          bValue1 -= aBase[-2*BLOCKING8+1]*bBase[1];
2242          bValue2 -= aBase[-3*BLOCKING8+1]*bBase[1];
2243          bValue3 -= aBase[-4*BLOCKING8+1]*bBase[1];
2244          bValue0 -= aBase[-1*BLOCKING8+0]*bBase[0];
2245          bValue1 -= aBase[-2*BLOCKING8+0]*bBase[0];
2246          bValue2 -= aBase[-3*BLOCKING8+0]*bBase[0];
2247          bValue3 -= aBase[-4*BLOCKING8+0]*bBase[0];
2248          bBase2[j]=bValue0;
2249          bBase2[j-1]=bValue1;
2250          bBase2[j-2]=bValue2;
2251          bBase2[j-3]=bValue3;
2252          aBase -= 4*BLOCKING8;
2253        }
2254      }
2255    }
2256  }
2257}
2258void CoinAbcDtrsm2(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2259{
2260  assert ((m&(BLOCKING8-1))==0);
2261  // 2 Left Lower Transpose Unit
2262  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2263    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2264    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2265      double * COIN_RESTRICT bBase = b+ii;
2266      double * COIN_RESTRICT bBase2 = bBase;
2267      const double * COIN_RESTRICT aBase=a+m*ii+(ii+BLOCKING8)*BLOCKING8;
2268      for (int i=BLOCKING8-1;i>=0;i--) {
2269        aBase -= BLOCKING8;
2270        for (int k=i+1;k<BLOCKING8;k++) {
2271          bBase2[i] -= aBase[k]*bBase2[k];
2272        }
2273      }
2274      for (int iii=ii-BLOCKING8;iii>=mm;iii-=BLOCKING8) {
2275        bBase2 -= BLOCKING8;
2276        assert (bBase2==b+iii);
2277        aBase -= m*BLOCKING8;
2278        const double * COIN_RESTRICT aBase2 = aBase+BLOCKING8X8;
2279        for (int i=BLOCKING8-1;i>=0;i--) {
2280          aBase2 -= BLOCKING8;
2281          for (int k=BLOCKING8-1;k>=0;k--) {
2282            bBase2[i] -= aBase2[k]*bBase[k];
2283          }
2284        }
2285      }
2286    }
2287    dtrsm2(kkk, 0, mm,  m, a, b);
2288  }
2289}
2290#define UNROLL_DTRSM3 16
2291#define CILK_DTRSM3 32
2292static void dtrsm3(int kkk, int first, int last,
2293                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2294{
2295  //int mm=CoinMin(kkk+UNROLL_DTRSM3*BLOCKING8,m);
2296  assert ((last-first)%BLOCKING8==0);
2297  if (last-first>CILK_DTRSM3) {
2298    int mid=((first+last)>>4)<<3;
2299    cilk_spawn dtrsm3(kkk,first,mid,m,a,b);
2300    dtrsm3(kkk,mid,last,m,a,b);
2301    cilk_sync;
2302  } else {
2303    for (int kk=0;kk<kkk;kk+=BLOCKING8) {
2304      for (int ii=first;ii<last;ii+=BLOCKING8) {
2305        double * COIN_RESTRICT bBase = b+ii;
2306        double * COIN_RESTRICT bBase2 = b+kk;
2307        const double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
2308        for (int j=0;j<BLOCKING8;j+=4) {
2309          double bValue0=bBase[j];
2310          double bValue1=bBase[j+1];
2311          double bValue2=bBase[j+2];
2312          double bValue3=bBase[j+3];
2313          bValue0 -= aBase[0]*bBase2[0];
2314          bValue1 -= aBase[1*BLOCKING8+0]*bBase2[0];
2315          bValue2 -= aBase[2*BLOCKING8+0]*bBase2[0];
2316          bValue3 -= aBase[3*BLOCKING8+0]*bBase2[0];
2317          bValue0 -= aBase[1]*bBase2[1];
2318          bValue1 -= aBase[1*BLOCKING8+1]*bBase2[1];
2319          bValue2 -= aBase[2*BLOCKING8+1]*bBase2[1];
2320          bValue3 -= aBase[3*BLOCKING8+1]*bBase2[1];
2321          bValue0 -= aBase[2]*bBase2[2];
2322          bValue1 -= aBase[1*BLOCKING8+2]*bBase2[2];
2323          bValue2 -= aBase[2*BLOCKING8+2]*bBase2[2];
2324          bValue3 -= aBase[3*BLOCKING8+2]*bBase2[2];
2325          bValue0 -= aBase[3]*bBase2[3];
2326          bValue1 -= aBase[1*BLOCKING8+3]*bBase2[3];
2327          bValue2 -= aBase[2*BLOCKING8+3]*bBase2[3];
2328          bValue3 -= aBase[3*BLOCKING8+3]*bBase2[3];
2329          bValue0 -= aBase[4]*bBase2[4];
2330          bValue1 -= aBase[1*BLOCKING8+4]*bBase2[4];
2331          bValue2 -= aBase[2*BLOCKING8+4]*bBase2[4];
2332          bValue3 -= aBase[3*BLOCKING8+4]*bBase2[4];
2333          bValue0 -= aBase[5]*bBase2[5];
2334          bValue1 -= aBase[1*BLOCKING8+5]*bBase2[5];
2335          bValue2 -= aBase[2*BLOCKING8+5]*bBase2[5];
2336          bValue3 -= aBase[3*BLOCKING8+5]*bBase2[5];
2337          bValue0 -= aBase[6]*bBase2[6];
2338          bValue1 -= aBase[1*BLOCKING8+6]*bBase2[6];
2339          bValue2 -= aBase[2*BLOCKING8+7]*bBase2[7];
2340          bValue3 -= aBase[3*BLOCKING8+6]*bBase2[6];
2341          bValue0 -= aBase[7]*bBase2[7];
2342          bValue1 -= aBase[1*BLOCKING8+7]*bBase2[7];
2343          bValue2 -= aBase[2*BLOCKING8+6]*bBase2[6];
2344          bValue3 -= aBase[3*BLOCKING8+7]*bBase2[7];
2345          bBase[j]=bValue0;
2346          bBase[j+1]=bValue1;
2347          bBase[j+2]=bValue2;
2348          bBase[j+3]=bValue3;
2349          aBase += 4*BLOCKING8;
2350        }
2351      }
2352    }
2353  }
2354}
2355void CoinAbcDtrsm3(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2356{
2357  assert ((m&(BLOCKING8-1))==0);
2358  // 3 Left Upper Transpose NonUnit
2359  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM3*BLOCKING8) {
2360    int mm=CoinMin(m,kkk+UNROLL_DTRSM3*BLOCKING8);
2361    dtrsm3(kkk, kkk, mm,  m, a, b);
2362    for (int ii=kkk;ii<mm;ii+=BLOCKING8) {
2363      double * COIN_RESTRICT bBase = b+ii;
2364      for (int kk=kkk;kk<ii;kk+=BLOCKING8) {
2365        double * COIN_RESTRICT bBase2 = b+kk;
2366        const double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
2367        for (int i=0;i<BLOCKING8;i++) {
2368          for (int k=0;k<BLOCKING8;k++) {
2369            bBase[i] -= aBase[k]*bBase2[k];
2370          }
2371          aBase += BLOCKING8;
2372        }
2373      }
2374      const double * COIN_RESTRICT aBase=a+ii*m+ii*BLOCKING8;
2375      for (int i=0;i<BLOCKING8;i++) {
2376        for (int k=0;k<i;k++) {
2377          bBase[i] -= aBase[k]*bBase[k];
2378        }
2379        bBase[i] = bBase[i]*aBase[i];
2380        aBase += BLOCKING8;
2381      }
2382    }
2383  }
2384}
2385static void 
2386CoinAbcDlaswp(int n, double * COIN_RESTRICT a,int lda,int start,int end, int * ipiv)
2387{
2388  assert ((n&(BLOCKING8-1))==0);
2389  double * COIN_RESTRICT aThis3 = a;
2390  for (int j=0;j<n;j+=BLOCKING8) {
2391    for (int i=start;i<end;i++) {
2392      int ip=ipiv[i];
2393      if (ip!=i) {
2394        double * COIN_RESTRICT aTop=aThis3+j*lda;
2395        int iBias = ip/BLOCKING8;
2396        ip-=iBias*BLOCKING8;
2397        double * COIN_RESTRICT aNotTop=aTop+iBias*BLOCKING8X8;
2398        iBias = i/BLOCKING8;
2399        int i2=i-iBias*BLOCKING8;
2400        aTop += iBias*BLOCKING8X8;
2401        for (int k=0;k<BLOCKING8;k++) {
2402          double temp=aTop[i2];
2403          aTop[i2]=aNotTop[ip];
2404          aNotTop[ip]=temp;
2405          aTop += BLOCKING8;
2406          aNotTop += BLOCKING8;
2407        }
2408      }
2409    }
2410  }
2411}
2412extern void CoinAbcDgemm(int m, int n, int k, double * COIN_RESTRICT a,int lda,
2413                          double * COIN_RESTRICT b,double * COIN_RESTRICT c
2414#if ABC_PARALLEL==2
2415                          ,int parallelMode
2416#endif
2417                         );
2418static int CoinAbcDgetrf2(int m, int n, double * COIN_RESTRICT a, int * ipiv)
2419{
2420  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
2421  assert (n==BLOCKING8);
2422  int dimension=CoinMin(m,n);
2423  /* entry for column j and row i (when multiple of BLOCKING8)
2424     is at aBlocked+j*m+i*BLOCKING8
2425  */
2426  assert (dimension==BLOCKING8);
2427  //printf("Dgetrf2 m=%d n=%d\n",m,n);
2428  double * COIN_RESTRICT aThis3 = a;
2429  for (int j=0;j<dimension;j++) {
2430    int pivotRow=-1;
2431    double largest=0.0;
2432    double * COIN_RESTRICT aThis2 = aThis3+j*BLOCKING8;
2433    // this block
2434    int pivotRow2=-1;
2435    for (int i=j;i<BLOCKING8;i++) {
2436      double value=fabs(aThis2[i]);
2437      if (value>largest) {
2438        largest=value;
2439        pivotRow2=i;
2440      }
2441    }
2442    // other blocks
2443    int iBias=0;
2444    for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
2445      aThis2+=BLOCKING8*BLOCKING8;
2446      for (int i=0;i<BLOCKING8;i++) {
2447        double value=fabs(aThis2[i]);
2448        if (value>largest) {
2449          largest=value;
2450          pivotRow2=i;
2451          iBias=ii;
2452        }
2453      }
2454    }
2455    pivotRow=pivotRow2+iBias;
2456    ipiv[j]=pivotRow;
2457    if (largest) {
2458      if (j!=pivotRow) {
2459        double * COIN_RESTRICT aTop=aThis3;
2460        double * COIN_RESTRICT aNotTop=aThis3+iBias*BLOCKING8;
2461        for (int i=0;i<n;i++) {
2462          double value = aTop[j];
2463          aTop[j]=aNotTop[pivotRow2];
2464          aNotTop[pivotRow2]=value;
2465          aTop += BLOCKING8;
2466          aNotTop += BLOCKING8;
2467        }
2468      }
2469      aThis2 = aThis3+j*BLOCKING8;
2470      double pivotMultiplier = 1.0 / aThis2[j];
2471      aThis2[j] = pivotMultiplier;
2472      // Do L
2473      // this block
2474      for (int i=j+1;i<BLOCKING8;i++) 
2475        aThis2[i] *= pivotMultiplier;
2476      // other blocks
2477      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
2478        aThis2+=BLOCKING8*BLOCKING8;
2479        for (int i=0;i<BLOCKING8;i++) 
2480          aThis2[i] *= pivotMultiplier;
2481      }
2482      // update rest
2483      double * COIN_RESTRICT aPut;
2484      aThis2 = aThis3+j*BLOCKING8;
2485      aPut = aThis2+BLOCKING8;
2486      double * COIN_RESTRICT aPut2 = aPut;
2487      // this block
2488      for (int i=j+1;i<BLOCKING8;i++) {
2489        double value = aPut[j];
2490        for (int k=j+1;k<BLOCKING8;k++) 
2491          aPut[k] -= value*aThis2[k];
2492        aPut += BLOCKING8;
2493      }
2494      // other blocks
2495      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
2496        aThis2 += BLOCKING8*BLOCKING8;
2497        aPut = aThis2+BLOCKING8;
2498        double * COIN_RESTRICT aPut2a = aPut2;
2499        for (int i=j+1;i<BLOCKING8;i++) {
2500          double value = aPut2a[j];
2501          for (int k=0;k<BLOCKING8;k++) 
2502            aPut[k] -= value*aThis2[k];
2503          aPut += BLOCKING8;
2504          aPut2a += BLOCKING8;
2505        }
2506      }
2507    } else {
2508      return 1;
2509    }
2510  }
2511  return 0;
2512}
2513
2514int 
2515CoinAbcDgetrf(int m, int n, double * COIN_RESTRICT a, int lda, int * ipiv
2516#if ABC_PARALLEL==2
2517                          ,int parallelMode
2518#endif
2519)
2520{
2521  assert (m==n);
2522  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
2523  if (m<BLOCKING8) {
2524    return CoinAbcDgetrf2(m,n,a,ipiv);
2525  } else {
2526    for (int j=0;j<n;j+=BLOCKING8) {
2527      int start=j;
2528      int newSize=CoinMin(BLOCKING8,n-j);
2529      int end=j+newSize;
2530      int returnCode=CoinAbcDgetrf2(m-start,newSize,a+(start*lda+start*BLOCKING8),
2531                                     ipiv+start);
2532      if (!returnCode) {
2533        // adjust
2534        for (int k=start;k<end;k++)
2535          ipiv[k]+=start;
2536        // swap 0<start
2537        CoinAbcDlaswp(start,a,lda,start,end,ipiv);
2538        if (end<n) {
2539          // swap >=end
2540          CoinAbcDlaswp(n-end,a+end*lda,lda,start,end,ipiv);
2541          CoinAbcDtrsmFactor(newSize,n-end,a+(start*lda+start*BLOCKING8),lda);
2542          CoinAbcDgemm(n-end,n-end,newSize,
2543                        a+start*lda+end*BLOCKING8,lda,
2544                        a+end*lda+start*BLOCKING8,a+end*lda+end*BLOCKING8
2545#if ABC_PARALLEL==2
2546                          ,parallelMode
2547#endif
2548                        );
2549        }
2550      } else {
2551        return returnCode;
2552      }
2553    }
2554  }
2555  return 0;
2556}
2557void
2558CoinAbcDgetrs(char trans,int m, double * COIN_RESTRICT a, double * COIN_RESTRICT work)
2559{
2560  assert ((m&(BLOCKING8-1))==0);
2561  if (trans=='N') {
2562    //CoinAbcDlaswp1(work,m,ipiv);
2563    CoinAbcDtrsm0(m,a,work);
2564    CoinAbcDtrsm1(m,a,work);
2565  } else {
2566    assert (trans=='T');
2567    CoinAbcDtrsm3(m,a,work);
2568    CoinAbcDtrsm2(m,a,work);
2569    //CoinAbcDlaswp1Backwards(work,m,ipiv);
2570  }
2571}
2572#ifdef ABC_LONG_FACTORIZATION
2573/// ****** Start long double version
2574/* type
2575   0 Left Lower NoTranspose Unit
2576   1 Left Upper NoTranspose NonUnit
2577   2 Left Lower Transpose Unit
2578   3 Left Upper Transpose NonUnit
2579*/
2580static void CoinAbcDtrsmFactor(int m, int n, long double * COIN_RESTRICT a,int lda)
2581{
2582  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
2583  assert (m==BLOCKING8);
2584  // 0 Left Lower NoTranspose Unit
2585  /* entry for column j and row i (when multiple of BLOCKING8)
2586     is at aBlocked+j*m+i*BLOCKING8
2587  */
2588  long double * COIN_RESTRICT aBase2 = a;
2589  long double * COIN_RESTRICT bBase2 = aBase2+lda*BLOCKING8;
2590  for (int jj=0;jj<n;jj+=BLOCKING8) {
2591    long double * COIN_RESTRICT bBase = bBase2;
2592    for (int j=jj;j<jj+BLOCKING8;j++) {
2593#if 0
2594      long double * COIN_RESTRICT aBase = aBase2;
2595      for (int k=0;k<BLOCKING8;k++) {
2596        long double bValue = bBase[k];
2597        if (bValue) {
2598          for (int i=k+1;i<BLOCKING8;i++) {
2599            bBase[i]-=bValue*aBase[i];
2600          }
2601        }
2602        aBase+=BLOCKING8;
2603      }
2604#else
2605      // unroll - stay in registers - don't test for zeros
2606      assert (BLOCKING8==8);
2607      long double bValue0 = bBase[0];
2608      long double bValue1 = bBase[1];
2609      long double bValue2 = bBase[2];
2610      long double bValue3 = bBase[3];
2611      long double bValue4 = bBase[4];
2612      long double bValue5 = bBase[5];
2613      long double bValue6 = bBase[6];
2614      long double bValue7 = bBase[7];
2615      bValue1-=bValue0*a[1+0*BLOCKING8];
2616      bBase[1] = bValue1;
2617      bValue2-=bValue0*a[2+0*BLOCKING8];
2618      bValue3-=bValue0*a[3+0*BLOCKING8];
2619      bValue4-=bValue0*a[4+0*BLOCKING8];
2620      bValue5-=bValue0*a[5+0*BLOCKING8];
2621      bValue6-=bValue0*a[6+0*BLOCKING8];
2622      bValue7-=bValue0*a[7+0*BLOCKING8];
2623      bValue2-=bValue1*a[2+1*BLOCKING8];
2624      bBase[2] = bValue2;
2625      bValue3-=bValue1*a[3+1*BLOCKING8];
2626      bValue4-=bValue1*a[4+1*BLOCKING8];
2627      bValue5-=bValue1*a[5+1*BLOCKING8];
2628      bValue6-=bValue1*a[6+1*BLOCKING8];
2629      bValue7-=bValue1*a[7+1*BLOCKING8];
2630      bValue3-=bValue2*a[3+2*BLOCKING8];
2631      bBase[3] = bValue3;
2632      bValue4-=bValue2*a[4+2*BLOCKING8];
2633      bValue5-=bValue2*a[5+2*BLOCKING8];
2634      bValue6-=bValue2*a[6+2*BLOCKING8];
2635      bValue7-=bValue2*a[7+2*BLOCKING8];
2636      bValue4-=bValue3*a[4+3*BLOCKING8];
2637      bBase[4] = bValue4;
2638      bValue5-=bValue3*a[5+3*BLOCKING8];
2639      bValue6-=bValue3*a[6+3*BLOCKING8];
2640      bValue7-=bValue3*a[7+3*BLOCKING8];
2641      bValue5-=bValue4*a[5+4*BLOCKING8];
2642      bBase[5] = bValue5;
2643      bValue6-=bValue4*a[6+4*BLOCKING8];
2644      bValue7-=bValue4*a[7+4*BLOCKING8];
2645      bValue6-=bValue5*a[6+5*BLOCKING8];
2646      bBase[6] = bValue6;
2647      bValue7-=bValue5*a[7+5*BLOCKING8];
2648      bValue7-=bValue6*a[7+6*BLOCKING8];
2649      bBase[7] = bValue7;
2650#endif
2651      bBase += BLOCKING8;
2652    }
2653    bBase2 +=lda*BLOCKING8; 
2654  }
2655}
2656#define UNROLL_DTRSM 16
2657#define CILK_DTRSM 32
2658static void dtrsm0(int kkk, int first, int last,
2659                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2660{
2661  int mm=CoinMin(kkk+UNROLL_DTRSM*BLOCKING8,m);
2662  assert ((last-first)%BLOCKING8==0);
2663  if (last-first>CILK_DTRSM) {
2664    int mid=((first+last)>>4)<<3;
2665    cilk_spawn dtrsm0(kkk,first,mid,m,a,b);
2666    dtrsm0(kkk,mid,last,m,a,b);
2667    cilk_sync;
2668  } else {
2669    const long double * COIN_RESTRICT aBaseA = a+UNROLL_DTRSM*BLOCKING8X8+kkk*BLOCKING8;
2670    aBaseA += (first-mm)*BLOCKING8-BLOCKING8X8;
2671    aBaseA += m*kkk;
2672    for (int ii=first;ii<last;ii+=BLOCKING8) {
2673      aBaseA += BLOCKING8X8;
2674      const long double * COIN_RESTRICT aBaseB=aBaseA;
2675      long double * COIN_RESTRICT bBaseI = b+ii;
2676      for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
2677        long double * COIN_RESTRICT bBase = b+kk;
2678        const long double * COIN_RESTRICT aBase2 = aBaseB;//a+UNROLL_DTRSM*BLOCKING8X8+m*kk+kkk*BLOCKING8;
2679        //aBase2 += (ii-mm)*BLOCKING8;
2680        //assert (aBase2==aBaseB);
2681        aBaseB+=m*BLOCKING8;
2682#if AVX2 !=2
2683#define ALTERNATE_INNER
2684#ifndef ALTERNATE_INNER
2685        for (int k=0;k<BLOCKING8;k++) {
2686          //coin_prefetch_const(aBase2+BLOCKING8);             
2687          for (int i=0;i<BLOCKING8;i++) {
2688            bBaseI[i] -= bBase[k]*aBase2[i];
2689          }
2690          aBase2 += BLOCKING8;
2691        }
2692#else
2693        long double b0=bBaseI[0];
2694        long double b1=bBaseI[1];
2695        long double b2=bBaseI[2];
2696        long double b3=bBaseI[3];
2697        long double b4=bBaseI[4];
2698        long double b5=bBaseI[5];
2699        long double b6=bBaseI[6];
2700        long double b7=bBaseI[7];
2701        for (int k=0;k<BLOCKING8;k++) {
2702          //coin_prefetch_const(aBase2+BLOCKING8);             
2703          long double bValue=bBase[k];
2704          b0 -= bValue * aBase2[0];
2705          b1 -= bValue * aBase2[1];
2706          b2 -= bValue * aBase2[2];
2707          b3 -= bValue * aBase2[3];
2708          b4 -= bValue * aBase2[4];
2709          b5 -= bValue * aBase2[5];
2710          b6 -= bValue * aBase2[6];
2711          b7 -= bValue * aBase2[7];
2712          aBase2 += BLOCKING8;
2713        }
2714        bBaseI[0]=b0;
2715        bBaseI[1]=b1;
2716        bBaseI[2]=b2;
2717        bBaseI[3]=b3;
2718        bBaseI[4]=b4;
2719        bBaseI[5]=b5;
2720        bBaseI[6]=b6;
2721        bBaseI[7]=b7;
2722#endif
2723#else
2724        __m256d b0=_mm256_load_pd(bBaseI);
2725        __m256d b1=_mm256_load_pd(bBaseI+4);
2726        for (int k=0;k<BLOCKING8;k++) {
2727          //coin_prefetch_const(aBase2+BLOCKING8);             
2728          __m256d bb = _mm256_broadcast_sd(bBase+k);
2729          __m256d a0 = _mm256_load_pd(aBase2);
2730          b0 -= a0*bb;
2731          __m256d a1 = _mm256_load_pd(aBase2+4);
2732          b1 -= a1*bb;
2733          aBase2 += BLOCKING8;
2734        }
2735        //_mm256_store_pd ((bBaseI), (b0));
2736        *reinterpret_cast<__m256d *>(bBaseI)=b0;
2737        //_mm256_store_pd (bBaseI+4, b1);
2738        *reinterpret_cast<__m256d *>(bBaseI+4)=b1;
2739#endif
2740      }
2741    }
2742  }
2743}
2744void CoinAbcDtrsm0(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2745{
2746  assert ((m&(BLOCKING8-1))==0);
2747  // 0 Left Lower NoTranspose Unit
2748  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM*BLOCKING8) {
2749    int mm=CoinMin(m,kkk+UNROLL_DTRSM*BLOCKING8);
2750    for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
2751      const long double * COIN_RESTRICT aBase2 = a+kk*(m+BLOCKING8);
2752      long double * COIN_RESTRICT bBase = b+kk;
2753      for (int k=0;k<BLOCKING8;k++) {
2754        for (int i=k+1;i<BLOCKING8;i++) {
2755          bBase[i] -= bBase[k]*aBase2[i];
2756        }
2757        aBase2 += BLOCKING8;
2758      }
2759      for (int ii=kk+BLOCKING8;ii<mm;ii+=BLOCKING8) {
2760        long double * COIN_RESTRICT bBaseI = b+ii;
2761        for (int k=0;k<BLOCKING8;k++) {
2762          //coin_prefetch_const(aBase2+BLOCKING8);             
2763          for (int i=0;i<BLOCKING8;i++) {
2764            bBaseI[i] -= bBase[k]*aBase2[i];
2765          }
2766          aBase2 += BLOCKING8;
2767        }
2768      }
2769    }
2770    dtrsm0(kkk,mm,m,m,a,b);
2771  }
2772}
2773static void dtrsm1(int kkk, int first, int last,
2774                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2775{
2776  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2777  assert ((last-first)%BLOCKING8==0);
2778  if (last-first>CILK_DTRSM) {
2779    int mid=((first+last)>>4)<<3;
2780    cilk_spawn dtrsm1(kkk,first,mid,m,a,b);
2781    dtrsm1(kkk,mid,last,m,a,b);
2782    cilk_sync;
2783  } else {
2784    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2785      long double * COIN_RESTRICT bBase2 = b+iii;
2786      const long double * COIN_RESTRICT aBaseA=
2787        a+BLOCKING8X8+BLOCKING8*iii;
2788      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2789        long double * COIN_RESTRICT bBase = b+ii;
2790        const long double * COIN_RESTRICT aBase=aBaseA+m*ii;
2791#if AVX2 !=2
2792#ifndef ALTERNATE_INNER
2793        for (int i=BLOCKING8-1;i>=0;i--) {
2794          aBase -= BLOCKING8;
2795          //coin_prefetch_const(aBase-BLOCKING8);               
2796          for (int k=BLOCKING8-1;k>=0;k--) {
2797            bBase2[k] -= bBase[i]*aBase[k];
2798          }
2799        }
2800#else
2801        long double b0=bBase2[0];
2802        long double b1=bBase2[1];
2803        long double b2=bBase2[2];
2804        long double b3=bBase2[3];
2805        long double b4=bBase2[4];
2806        long double b5=bBase2[5];
2807        long double b6=bBase2[6];
2808        long double b7=bBase2[7];
2809        for (int i=BLOCKING8-1;i>=0;i--) {
2810          aBase -= BLOCKING8;
2811          //coin_prefetch_const(aBase-BLOCKING8);               
2812          long double bValue=bBase[i];
2813          b0 -= bValue * aBase[0];
2814          b1 -= bValue * aBase[1];
2815          b2 -= bValue * aBase[2];
2816          b3 -= bValue * aBase[3];
2817          b4 -= bValue * aBase[4];
2818          b5 -= bValue * aBase[5];
2819          b6 -= bValue * aBase[6];
2820          b7 -= bValue * aBase[7];
2821        }
2822        bBase2[0]=b0;
2823        bBase2[1]=b1;
2824        bBase2[2]=b2;
2825        bBase2[3]=b3;
2826        bBase2[4]=b4;
2827        bBase2[5]=b5;
2828        bBase2[6]=b6;
2829        bBase2[7]=b7;
2830#endif
2831#else
2832        __m256d b0=_mm256_load_pd(bBase2);
2833        __m256d b1=_mm256_load_pd(bBase2+4);
2834        for (int i=BLOCKING8-1;i>=0;i--) {
2835          aBase -= BLOCKING8;
2836          //coin_prefetch_const(aBase5-BLOCKING8);             
2837          __m256d bb = _mm256_broadcast_sd(bBase+i);
2838          __m256d a0 = _mm256_load_pd(aBase);
2839          b0 -= a0*bb;
2840          __m256d a1 = _mm256_load_pd(aBase+4);
2841          b1 -= a1*bb;
2842        }
2843        //_mm256_store_pd (bBase2, b0);
2844        *reinterpret_cast<__m256d *>(bBase2)=b0;
2845        //_mm256_store_pd (bBase2+4, b1);
2846        *reinterpret_cast<__m256d *>(bBase2+4)=b1;
2847#endif
2848      }
2849    }
2850  }
2851}
2852void CoinAbcDtrsm1(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2853{
2854  assert ((m&(BLOCKING8-1))==0);
2855  // 1 Left Upper NoTranspose NonUnit
2856  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2857    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2858    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2859      const long double * COIN_RESTRICT aBase = a+m*ii+ii*BLOCKING8;
2860      long double * COIN_RESTRICT bBase = b+ii;
2861      for (int i=BLOCKING8-1;i>=0;i--) {
2862        bBase[i] *= aBase[i*(BLOCKING8+1)];
2863        for (int k=i-1;k>=0;k--) {
2864          bBase[k] -= bBase[i]*aBase[k+i*BLOCKING8];
2865        }
2866      }
2867      long double * COIN_RESTRICT bBase2 = bBase;
2868      for (int iii=ii;iii>mm;iii-=BLOCKING8) {
2869        bBase2 -= BLOCKING8;
2870        for (int i=BLOCKING8-1;i>=0;i--) {
2871          aBase -= BLOCKING8;
2872          coin_prefetch_const(aBase-BLOCKING8);         
2873          for (int k=BLOCKING8-1;k>=0;k--) {
2874            bBase2[k] -= bBase[i]*aBase[k];
2875          }
2876        }
2877      }
2878    }
2879    dtrsm1(kkk, 0, mm,  m, a, b);
2880  }
2881}
2882static void dtrsm2(int kkk, int first, int last,
2883                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2884{
2885  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2886  assert ((last-first)%BLOCKING8==0);
2887  if (last-first>CILK_DTRSM) {
2888    int mid=((first+last)>>4)<<3;
2889    cilk_spawn dtrsm2(kkk,first,mid,m,a,b);
2890    dtrsm2(kkk,mid,last,m,a,b);
2891    cilk_sync;
2892  } else {
2893    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2894      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2895        const long double * COIN_RESTRICT bBase = b+ii;
2896        long double * COIN_RESTRICT bBase2=b+iii;
2897        const long double * COIN_RESTRICT aBase=a+ii*BLOCKING8+m*iii+BLOCKING8X8;
2898        for (int j=BLOCKING8-1;j>=0;j-=4) {
2899          long double bValue0=bBase2[j];
2900          long double bValue1=bBase2[j-1];
2901          long double bValue2=bBase2[j-2];
2902          long double bValue3=bBase2[j-3];
2903          bValue0 -= aBase[-1*BLOCKING8+7]*bBase[7];
2904          bValue1 -= aBase[-2*BLOCKING8+7]*bBase[7];
2905          bValue2 -= aBase[-3*BLOCKING8+7]*bBase[7];
2906          bValue3 -= aBase[-4*BLOCKING8+7]*bBase[7];
2907          bValue0 -= aBase[-1*BLOCKING8+6]*bBase[6];
2908          bValue1 -= aBase[-2*BLOCKING8+6]*bBase[6];
2909          bValue2 -= aBase[-3*BLOCKING8+6]*bBase[6];
2910          bValue3 -= aBase[-4*BLOCKING8+6]*bBase[6];
2911          bValue0 -= aBase[-1*BLOCKING8+5]*bBase[5];
2912          bValue1 -= aBase[-2*BLOCKING8+5]*bBase[5];
2913          bValue2 -= aBase[-3*BLOCKING8+5]*bBase[5];
2914          bValue3 -= aBase[-4*BLOCKING8+5]*bBase[5];
2915          bValue0 -= aBase[-1*BLOCKING8+4]*bBase[4];
2916          bValue1 -= aBase[-2*BLOCKING8+4]*bBase[4];
2917          bValue2 -= aBase[-3*BLOCKING8+4]*bBase[4];
2918          bValue3 -= aBase[-4*BLOCKING8+4]*bBase[4];
2919          bValue0 -= aBase[-1*BLOCKING8+3]*bBase[3];
2920          bValue1 -= aBase[-2*BLOCKING8+3]*bBase[3];
2921          bValue2 -= aBase[-3*BLOCKING8+3]*bBase[3];
2922          bValue3 -= aBase[-4*BLOCKING8+3]*bBase[3];
2923          bValue0 -= aBase[-1*BLOCKING8+2]*bBase[2];
2924          bValue1 -= aBase[-2*BLOCKING8+2]*bBase[2];
2925          bValue2 -= aBase[-3*BLOCKING8+2]*bBase[2];
2926          bValue3 -= aBase[-4*BLOCKING8+2]*bBase[2];
2927          bValue0 -= aBase[-1*BLOCKING8+1]*bBase[1];
2928          bValue1 -= aBase[-2*BLOCKING8+1]*bBase[1];
2929          bValue2 -= aBase[-3*BLOCKING8+1]*bBase[1];
2930          bValue3 -= aBase[-4*BLOCKING8+1]*bBase[1];
2931          bValue0 -= aBase[-1*BLOCKING8+0]*bBase[0];
2932          bValue1 -= aBase[-2*BLOCKING8+0]*bBase[0];
2933          bValue2 -= aBase[-3*BLOCKING8+0]*bBase[0];
2934          bValue3 -= aBase[-4*BLOCKING8+0]*bBase[0];
2935          bBase2[j]=bValue0;
2936          bBase2[j-1]=bValue1;
2937          bBase2[j-2]=bValue2;
2938          bBase2[j-3]=bValue3;
2939          aBase -= 4*BLOCKING8;
2940        }
2941      }
2942    }
2943  }
2944}
2945void CoinAbcDtrsm2(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2946{
2947  assert ((m&(BLOCKING8-1))==0);
2948  // 2 Left Lower Transpose Unit
2949  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2950    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2951    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2952      long double * COIN_RESTRICT bBase = b+ii;
2953      long double * COIN_RESTRICT bBase2 = bBase;
2954      const long double * COIN_RESTRICT aBase=a+m*ii+(ii+BLOCKING8)*BLOCKING8;
2955      for (int i=BLOCKING8-1;i>=0;i--) {
2956        aBase -= BLOCKING8;
2957        for (int k=i+1;k<BLOCKING8;k++) {
2958          bBase2[i] -= aBase[k]*bBase2[k];
2959        }
2960      }
2961      for (int iii=ii-BLOCKING8;iii>=mm;iii-=BLOCKING8) {
2962        bBase2 -= BLOCKING8;
2963        assert (bBase2==b+iii);
2964        aBase -= m*BLOCKING8;
2965        const long double * COIN_RESTRICT aBase2 = aBase+BLOCKING8X8;
2966        for (int i=BLOCKING8-1;i>=0;i--) {
2967          aBase2 -= BLOCKING8;
2968          for (int k=BLOCKING8-1;k>=0;k--) {
2969            bBase2[i] -= aBase2[k]*bBase[k];
2970          }
2971        }
2972      }
2973    }
2974    dtrsm2(kkk, 0, mm,  m, a, b);
2975  }
2976}
2977#define UNROLL_DTRSM3 16
2978#define CILK_DTRSM3 32
2979static void dtrsm3(int kkk, int first, int last,
2980                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2981{
2982  //int mm=CoinMin(kkk+UNROLL_DTRSM3*BLOCKING8,m);
2983  assert ((last-first)%BLOCKING8==0);
2984  if (last-first>CILK_DTRSM3) {
2985    int mid=((first+last)>>4)<<3;
2986    cilk_spawn dtrsm3(kkk,first,mid,m,a,b);
2987    dtrsm3(kkk,mid,last,m,a,b);
2988    cilk_sync;
2989  } else {
2990    for (int kk=0;kk<kkk;kk+=BLOCKING8) {
2991      for (int ii=first;ii<last;ii+=BLOCKING8) {
2992        long double * COIN_RESTRICT bBase = b+ii;
2993        long double * COIN_RESTRICT bBase2 = b+kk;
2994        const long double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
2995        for (int j=0;j<BLOCKING8;j+=4) {
2996          long double bValue0=bBase[j];
2997          long double bValue1=bBase[j+1];
2998          long double bValue2=bBase[j+2];
2999          long double bValue3=bBase[j+3];
3000          bValue0 -= aBase[0]*bBase2[0];
3001          bValue1 -= aBase[1*BLOCKING8+0]*bBase2[0];
3002          bValue2 -= aBase[2*BLOCKING8+0]*bBase2[0];
3003          bValue3 -= aBase[3*BLOCKING8+0]*bBase2[0];
3004          bValue0 -= aBase[1]*bBase2[1];
3005          bValue1 -= aBase[1*BLOCKING8+1]*bBase2[1];
3006          bValue2 -= aBase[2*BLOCKING8+1]*bBase2[1];
3007          bValue3 -= aBase[3*BLOCKING8+1]*bBase2[1];
3008          bValue0 -= aBase[2]*bBase2[2];
3009          bValue1 -= aBase[1*BLOCKING8+2]*bBase2[2];
3010          bValue2 -= aBase[2*BLOCKING8+2]*bBase2[2];
3011          bValue3 -= aBase[3*BLOCKING8+2]*bBase2[2];
3012          bValue0 -= aBase[3]*bBase2[3];
3013          bValue1 -= aBase[1*BLOCKING8+3]*bBase2[3];
3014          bValue2 -= aBase[2*BLOCKING8+3]*bBase2[3];
3015          bValue3 -= aBase[3*BLOCKING8+3]*bBase2[3];
3016          bValue0 -= aBase[4]*bBase2[4];
3017          bValue1 -= aBase[1*BLOCKING8+4]*bBase2[4];
3018          bValue2 -= aBase[2*BLOCKING8+4]*bBase2[4];
3019          bValue3 -= aBase[3*BLOCKING8+4]*bBase2[4];
3020          bValue0 -= aBase[5]*bBase2[5];
3021          bValue1 -= aBase[1*BLOCKING8+5]*bBase2[5];
3022          bValue2 -= aBase[2*BLOCKING8+5]*bBase2[5];
3023          bValue3 -= aBase[3*BLOCKING8+5]*bBase2[5];
3024          bValue0 -= aBase[6]*bBase2[6];
3025          bValue1 -= aBase[1*BLOCKING8+6]*bBase2[6];
3026          bValue2 -= aBase[2*BLOCKING8+7]*bBase2[7];
3027          bValue3 -= aBase[3*BLOCKING8+6]*bBase2[6];
3028          bValue0 -= aBase[7]*bBase2[7];
3029          bValue1 -= aBase[1*BLOCKING8+7]*bBase2[7];
3030          bValue2 -= aBase[2*BLOCKING8+6]*bBase2[6];
3031          bValue3 -= aBase[3*BLOCKING8+7]*bBase2[7];
3032          bBase[j]=bValue0;
3033          bBase[j+1]=bValue1;
3034          bBase[j+2]=bValue2;
3035          bBase[j+3]=bValue3;
3036          aBase += 4*BLOCKING8;
3037        }
3038      }
3039    }
3040  }
3041}
3042void CoinAbcDtrsm3(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
3043{
3044  assert ((m&(BLOCKING8-1))==0);
3045  // 3 Left Upper Transpose NonUnit
3046  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM3*BLOCKING8) {
3047    int mm=CoinMin(m,kkk+UNROLL_DTRSM3*BLOCKING8);
3048    dtrsm3(kkk, kkk, mm,  m, a, b);
3049    for (int ii=kkk;ii<mm;ii+=BLOCKING8) {
3050      long double * COIN_RESTRICT bBase = b+ii;
3051      for (int kk=kkk;kk<ii;kk+=BLOCKING8) {
3052        long double * COIN_RESTRICT bBase2 = b+kk;
3053        const long double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
3054        for (int i=0;i<BLOCKING8;i++) {
3055          for (int k=0;k<BLOCKING8;k++) {
3056            bBase[i] -= aBase[k]*bBase2[k];
3057          }
3058          aBase += BLOCKING8;
3059        }
3060      }
3061      const long double * COIN_RESTRICT aBase=a+ii*m+ii*BLOCKING8;
3062      for (int i=0;i<BLOCKING8;i++) {
3063        for (int k=0;k<i;k++) {
3064          bBase[i] -= aBase[k]*bBase[k];
3065        }
3066        bBase[i] = bBase[i]*aBase[i];
3067        aBase += BLOCKING8;
3068      }
3069    }
3070  }
3071}
3072static void 
3073CoinAbcDlaswp(int n, long double * COIN_RESTRICT a,int lda,int start,int end, int * ipiv)
3074{
3075  assert ((n&(BLOCKING8-1))==0);
3076  long double * COIN_RESTRICT aThis3 = a;
3077  for (int j=0;j<n;j+=BLOCKING8) {
3078    for (int i=start;i<end;i++) {
3079      int ip=ipiv[i];
3080      if (ip!=i) {
3081        long double * COIN_RESTRICT aTop=aThis3+j*lda;
3082        int iBias = ip/BLOCKING8;
3083        ip-=iBias*BLOCKING8;
3084        long double * COIN_RESTRICT aNotTop=aTop+iBias*BLOCKING8X8;
3085        iBias = i/BLOCKING8;
3086        int i2=i-iBias*BLOCKING8;
3087        aTop += iBias*BLOCKING8X8;
3088        for (int k=0;k<BLOCKING8;k++) {
3089          long double temp=aTop[i2];
3090          aTop[i2]=aNotTop[ip];
3091          aNotTop[ip]=temp;
3092          aTop += BLOCKING8;
3093          aNotTop += BLOCKING8;
3094        }
3095      }
3096    }
3097  }
3098}
3099extern void CoinAbcDgemm(int m, int n, int k, long double * COIN_RESTRICT a,int lda,
3100                          long double * COIN_RESTRICT b,long double * COIN_RESTRICT c
3101#if ABC_PARALLEL==2
3102                          ,int parallelMode
3103#endif
3104                         );
3105static int CoinAbcDgetrf2(int m, int n, long double * COIN_RESTRICT a, int * ipiv)
3106{
3107  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
3108  assert (n==BLOCKING8);
3109  int dimension=CoinMin(m,n);
3110  /* entry for column j and row i (when multiple of BLOCKING8)
3111     is at aBlocked+j*m+i*BLOCKING8
3112  */
3113  assert (dimension==BLOCKING8);
3114  //printf("Dgetrf2 m=%d n=%d\n",m,n);
3115  long double * COIN_RESTRICT aThis3 = a;
3116  for (int j=0;j<dimension;j++) {
3117    int pivotRow=-1;
3118    long double largest=0.0;
3119    long double * COIN_RESTRICT aThis2 = aThis3+j*BLOCKING8;
3120    // this block
3121    int pivotRow2=-1;
3122    for (int i=j;i<BLOCKING8;i++) {
3123      long double value=fabsl(aThis2[i]);
3124      if (value>largest) {
3125        largest=value;
3126        pivotRow2=i;
3127      }
3128    }
3129    // other blocks
3130    int iBias=0;
3131    for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
3132      aThis2+=BLOCKING8*BLOCKING8;
3133      for (int i=0;i<BLOCKING8;i++) {
3134        long double value=fabsl(aThis2[i]);
3135        if (value>largest) {
3136          largest=value;
3137          pivotRow2=i;
3138          iBias=ii;
3139        }
3140      }
3141    }
3142    pivotRow=pivotRow2+iBias;
3143    ipiv[j]=pivotRow;
3144    if (largest) {
3145      if (j!=pivotRow) {
3146        long double * COIN_RESTRICT aTop=aThis3;
3147        long double * COIN_RESTRICT aNotTop=aThis3+iBias*BLOCKING8;
3148        for (int i=0;i<n;i++) {
3149          long double value = aTop[j];
3150          aTop[j]=aNotTop[pivotRow2];
3151          aNotTop[pivotRow2]=value;
3152          aTop += BLOCKING8;
3153          aNotTop += BLOCKING8;
3154        }
3155      }
3156      aThis2 = aThis3+j*BLOCKING8;
3157      long double pivotMultiplier = 1.0 / aThis2[j];
3158      aThis2[j] = pivotMultiplier;
3159      // Do L
3160      // this block
3161      for (int i=j+1;i<BLOCKING8;i++) 
3162        aThis2[i] *= pivotMultiplier;
3163      // other blocks
3164      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
3165        aThis2+=BLOCKING8*BLOCKING8;
3166        for (int i=0;i<BLOCKING8;i++) 
3167          aThis2[i] *= pivotMultiplier;
3168      }
3169      // update rest
3170      long double * COIN_RESTRICT aPut;
3171      aThis2 = aThis3+j*BLOCKING8;
3172      aPut = aThis2+BLOCKING8;
3173      long double * COIN_RESTRICT aPut2 = aPut;
3174      // this block
3175      for (int i=j+1;i<BLOCKING8;i++) {
3176        long double value = aPut[j];
3177        for (int k=j+1;k<BLOCKING8;k++) 
3178          aPut[k] -= value*aThis2[k];
3179        aPut += BLOCKING8;
3180      }
3181      // other blocks
3182      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
3183        aThis2 += BLOCKING8*BLOCKING8;
3184        aPut = aThis2+BLOCKING8;
3185        long double * COIN_RESTRICT aPut2a = aPut2;
3186        for (int i=j+1;i<BLOCKING8;i++) {
3187          long double value = aPut2a[j];
3188          for (int k=0;k<BLOCKING8;k++) 
3189            aPut[k] -= value*aThis2[k];
3190          aPut += BLOCKING8;
3191          aPut2a += BLOCKING8;
3192        }
3193      }
3194    } else {
3195      return 1;
3196    }
3197  }
3198  return 0;
3199}
3200
3201int 
3202CoinAbcDgetrf(int m, int n, long double * COIN_RESTRICT a, int lda, int * ipiv
3203#if ABC_PARALLEL==2
3204                          ,int parallelMode
3205#endif
3206)
3207{
3208  assert (m==n);
3209  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
3210  if (m<BLOCKING8) {
3211    return CoinAbcDgetrf2(m,n,a,ipiv);
3212  } else {
3213    for (int j=0;j<n;j+=BLOCKING8) {
3214      int start=j;
3215      int newSize=CoinMin(BLOCKING8,n-j);
3216      int end=j+newSize;
3217      int returnCode=CoinAbcDgetrf2(m-start,newSize,a+(start*lda+start*BLOCKING8),
3218                                     ipiv+start);
3219      if (!returnCode) {
3220        // adjust
3221        for (int k=start;k<end;k++)
3222          ipiv[k]+=start;
3223        // swap 0<start
3224        CoinAbcDlaswp(start,a,lda,start,end,ipiv);
3225        if (end<n) {
3226          // swap >=end
3227          CoinAbcDlaswp(n-end,a+end*lda,lda,start,end,ipiv);
3228          CoinAbcDtrsmFactor(newSize,n-end,a+(start*lda+start*BLOCKING8),lda);
3229          CoinAbcDgemm(n-end,n-end,newSize,
3230                        a+start*lda+end*BLOCKING8,lda,
3231                        a+end*lda+start*BLOCKING8,a+end*lda+end*BLOCKING8
3232#if ABC_PARALLEL==2
3233                          ,parallelMode
3234#endif
3235                        );
3236        }
3237      } else {
3238        return returnCode;
3239      }
3240    }
3241  }
3242  return 0;
3243}
3244void
3245CoinAbcDgetrs(char trans,int m, long double * COIN_RESTRICT a, long double * COIN_RESTRICT work)
3246{
3247  assert ((m&(BLOCKING8-1))==0);
3248  if (trans=='N') {
3249    //CoinAbcDlaswp1(work,m,ipiv);
3250    CoinAbcDtrsm0(m,a,work);
3251    CoinAbcDtrsm1(m,a,work);
3252  } else {
3253    assert (trans=='T');
3254    CoinAbcDtrsm3(m,a,work);
3255    CoinAbcDtrsm2(m,a,work);
3256    //CoinAbcDlaswp1Backwards(work,m,ipiv);
3257  }
3258}
3259#endif
3260#endif
Note: See TracBrowser for help on using the repository browser.