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

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

fix a few problems with Aboca

File size: 96.5 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#if MOVE_REPLACE_PART1A > 0
1135  if (!abcFactorization_->usingFT()) {
1136#endif
1137    usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
1138    abcFactorization_->updateColumnTranspose(usefulArray_[arrayForBtran_]);
1139#if MOVE_REPLACE_PART1A > 0
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#if MOVE_REPLACE_PART1A <= 0
1181  cilk_spawn getTableauColumnPart2();
1182#if MOVE_REPLACE_PART1A == 0
1183  cilk_spawn checkReplacePart1();
1184#endif
1185  numberFlipped=flipBounds();
1186  cilk_sync;
1187#else
1188  if (abcFactorization_->usingFT()) {
1189    cilk_spawn getTableauColumnPart2();
1190    ftAlpha_ = cilk_spawn abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
1191    numberFlipped=flipBounds();
1192    cilk_sync;
1193  } else {
1194    cilk_spawn getTableauColumnPart2();
1195    numberFlipped=flipBounds();
1196    cilk_sync;
1197  }
1198#endif
1199  //usefulArray_[arrayForTableauRow_].compact();
1200  // returns 3 if skip this iteration and re-factorize
1201  /*
1202    return code
1203    0 - OK
1204    2 - flag something and skip
1205    3 - break and re-factorize
1206    4 - break and say infeasible
1207 */
1208  int returnCode=0;
1209  // amount primal will move
1210  movement_ = -dualOut_ * directionOut_ / alpha_;
1211  // see if update stable
1212#if ABC_NORMAL_DEBUG>3
1213  if ((handler_->logLevel() & 32)) {
1214    double ft=ftAlpha_*abcFactorization_->pivotRegion()[pivotRow_];
1215    double diff1=fabs(alpha_-btranAlpha_);
1216    double diff2=fabs(fabs(alpha_)-fabs(ft));
1217    double diff3=fabs(fabs(ft)-fabs(btranAlpha_));
1218    double largest=CoinMax(CoinMax(diff1,diff2),diff3);
1219    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n", 
1220           btranAlpha_, alpha_,ft,largest);
1221    if (largest>0.001*fabs(alpha_)) {
1222      printf("bad\n");
1223    }
1224  }
1225#endif
1226  int numberPivots=abcFactorization_->pivots();
1227  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
1228  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
1229  // if can't trust much and long way from optimal then relax
1230  if (largestPrimalError_ > 10.0)
1231    checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
1232  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1233      fabs(btranAlpha_ - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
1234    handler_->message(CLP_DUAL_CHECK, messages_)
1235      << btranAlpha_
1236      << alpha_
1237      << CoinMessageEol;
1238    // see with more relaxed criterion
1239    double test;
1240    if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
1241      test = 1.0e-1 * fabs(alpha_);
1242    else
1243      test = CoinMin(10.0*checkValue,1.0e-4 * (1.0 + fabs(alpha_)));
1244    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1245                     fabs(btranAlpha_ - alpha_) > test);
1246    double derror = CoinMin(fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),1.0)*0.9999e7;
1247    int error=0;
1248    while (derror>1.0) {
1249      error++;
1250      derror *= 0.1;
1251    }
1252    int frequency[8]={99999999,100,10,2,1,1,1,1};
1253    int newFactorFrequency=CoinMin(forceFactorization_,frequency[error]);
1254#if ABC_NORMAL_DEBUG>0
1255    if (newFactorFrequency<forceFactorization_)
1256      printf("Error of %g after %d pivots old forceFact %d now %d\n",
1257             fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_)),numberPivots,
1258             forceFactorization_,newFactorFrequency);
1259#endif
1260    if (!numberPivots&&fabs(btranAlpha_ - alpha_)/(1.0 + fabs(alpha_))>0.5e-6
1261        &&abcFactorization_->pivotTolerance()<0.5)
1262      abcFactorization_->saferTolerances (1.0, -1.03);
1263    forceFactorization_=newFactorFrequency;
1264
1265   
1266#if ABC_NORMAL_DEBUG>0
1267    if (fabs(btranAlpha_ + alpha_) < 1.0e-5*(1.0 + fabs(alpha_))) {
1268      printf("diff (%g,%g) %g check %g\n",btranAlpha_,alpha_,fabs(btranAlpha_-alpha_),checkValue*(1.0+fabs(alpha_)));
1269    }
1270#endif
1271    if (numberPivots) {
1272      if (needFlag&&numberPivots<10) {
1273        // need to reject something
1274        assert (sequenceOut_>=0);
1275        char x = isColumn(sequenceOut_) ? 'C' : 'R';
1276        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1277          << x << sequenceWithin(sequenceOut_)
1278          << CoinMessageEol;
1279#if ABC_NORMAL_DEBUG>0
1280        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
1281               btranAlpha_, alpha_,numberPivots);
1282#endif
1283        // Make safer?
1284        abcFactorization_->saferTolerances (1.0, -1.03);
1285        setFlagged(sequenceOut_);
1286        // so can't happen immediately again
1287        sequenceOut_=-1;
1288        lastBadIteration_ = numberIterations_; // say be more cautious
1289      }
1290      // Make safer?
1291      //if (numberPivots<5)
1292      //abcFactorization_->saferTolerances (-0.99, -1.01);
1293      problemStatus_ = -2; // factorize now
1294      returnCode = -2;
1295      stateOfIteration_=2;
1296    } else {
1297      if (needFlag) {
1298        assert (sequenceOut_>=0);
1299        // need to reject something
1300        char x = isColumn(sequenceOut_) ? 'C' : 'R';
1301        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1302          << x << sequenceWithin(sequenceOut_)
1303          << CoinMessageEol;
1304#if ABC_NORMAL_DEBUG>3
1305        printf("flag a %g %g\n", btranAlpha_, alpha_);
1306#endif
1307        setFlagged(sequenceOut_);
1308        // so can't happen immediately again
1309        sequenceOut_=-1;
1310        //abcProgress_.clearBadTimes();
1311        lastBadIteration_ = numberIterations_; // say be more cautious
1312        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
1313          //printf("I think should declare infeasible\n");
1314          problemStatus_ = 1;
1315          returnCode = 1;
1316          stateOfIteration_=2;
1317        } else {
1318          stateOfIteration_=1;
1319        }
1320        // Make safer?
1321        if (abcFactorization_->pivotTolerance()<0.999&&stateOfIteration_==1) {
1322          // change tolerance and re-invert
1323          abcFactorization_->saferTolerances (1.0, -1.03);
1324          problemStatus_ = -6; // factorize now
1325          returnCode = -2;
1326          stateOfIteration_=2;
1327        }
1328      }
1329    }
1330  }
1331  if (!stateOfIteration_) {
1332    // check update
1333    //int updateStatus =
1334/*
1335  returns
1336  0 - OK
1337  1 - take iteration then re-factorize
1338  2 - flag something and skip
1339  3 - break and re-factorize
1340  5 - take iteration then re-factorize because of memory
1341 */
1342    int status=checkReplace();
1343    if (status&&!returnCode)
1344      returnCode=status;
1345  }
1346  //clearArrays(arrayForFlipRhs_);
1347  if (stateOfIteration_&&numberFlipped) {
1348    //usefulArray_[arrayForTableauRow_].compact();
1349    // move variables back
1350    flipBack(numberFlipped);
1351  }
1352  // could choose average of all three
1353  return returnCode;
1354}
1355int 
1356AbcSimplexDual::whileIteratingParallel(int numberIterations)
1357{
1358  int returnCode=-1;
1359#ifdef EARLY_FACTORIZE
1360  int savePivot=-1;
1361  CoinIndexedVector & early = usefulArray_[ABC_NUMBER_USEFUL-1];
1362  int * pivotIndices = early.getIndices();
1363  double * pivotValues = early.denseVector();
1364  int pivotCountPosition=early.capacity()-1;
1365  int numberSaved=0;
1366  if (numberIterations==123456789)
1367    savePivot=pivotCountPosition;;
1368#endif
1369  numberIterations += numberIterations_;
1370  do {
1371#if ABC_DEBUG
1372    checkArrays();
1373#endif
1374    /*
1375      -1 - just move dual in values pass
1376      0 - take iteration
1377      1 - don't take but continue
1378      2 - don't take and break
1379    */
1380    stateOfIteration_=0;
1381    returnCode=-1;
1382    // put row of tableau in usefulArray[arrayForTableauRow_]
1383    /*
1384      Could
1385      a) copy [arrayForBtran_] and start off updateWeights earlier
1386      b) break transposeTimes into two and do after slack part
1387      c) also think if cleaner to go all way with update but just don't do final part
1388    */
1389    //upperTheta=COIN_DBL_MAX;
1390    double saveAcceptable=acceptablePivot_;
1391    if (largestPrimalError_>1.0e-5||largestDualError_>1.0e-5) {
1392    //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) {
1393      if (!abcFactorization_->pivots())
1394        acceptablePivot_*=1.0e2;
1395      else if (abcFactorization_->pivots()<5)
1396        acceptablePivot_*=1.0e1;
1397    }
1398#ifdef MOVE_UPDATE_WEIGHTS
1399    // copy btran across
1400    usefulArray_[5].copy(usefulArray_[arrayForBtran_]);
1401    cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[5]);;
1402#endif
1403    dualColumn1();
1404    acceptablePivot_=saveAcceptable;
1405    if ((stateOfProblem_&VALUES_PASS)!=0) {
1406      // see if can just move dual
1407      if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
1408        stateOfIteration_=-1;
1409      }
1410    }
1411    if (!stateOfIteration_) {
1412#ifndef MOVE_UPDATE_WEIGHTS
1413      cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[arrayForBtran_]);;
1414#endif
1415      // get sequenceIn_
1416      dualPivotColumn();
1417      //stateOfIteration_=0;
1418      if (sequenceIn_ >= 0) {
1419        // normal iteration
1420        // update the incoming column
1421        //arrayForReplaceColumn_=getAvailableArray();
1422        getTableauColumnFlipAndStartReplaceCilk();
1423        //usleep(1000);
1424      } else if (stateOfIteration_!=-1) {
1425        stateOfIteration_=2;
1426      }
1427    }
1428    cilk_sync;
1429    // Check event
1430    {
1431      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1432      if (status >= 0) {
1433        problemStatus_ = 5;
1434        secondaryStatus_ = ClpEventHandler::endOfIteration;
1435        returnCode = 4;
1436        stateOfIteration_=2;
1437      }
1438    }
1439    if (!stateOfIteration_) {
1440      //        assert (stateOfIteration_!=1);
1441      int whatNext=0;
1442      whatNext = housekeeping();
1443#ifdef EARLY_FACTORIZE
1444      if (savePivot>=0) {
1445        // save pivot
1446        pivotIndices[--savePivot]=sequenceOut_;
1447        pivotValues[savePivot]=alpha_;
1448        pivotIndices[--savePivot]=sequenceIn_;
1449        pivotValues[savePivot]=btranAlpha_;
1450        numberSaved++;
1451      }
1452#endif
1453      if (whatNext == 1) {
1454        problemStatus_ = -2; // refactorize
1455        usefulArray_[arrayForTableauRow_].compact();
1456      } else if (whatNext == 2) {
1457        // maximum iterations or equivalent
1458        problemStatus_ = 3;
1459        returnCode = 3;
1460        stateOfIteration_=2;
1461      }
1462#ifdef EARLY_FACTORIZE
1463      if (savePivot>=0) {
1464        // tell worker can update this
1465        pivotIndices[pivotCountPosition]=numberSaved;
1466      }
1467#endif
1468    } else {
1469      usefulArray_[arrayForTableauRow_].compact();
1470       if (stateOfIteration_<0) {
1471        // move dual in dual values pass
1472        theta_=abcDj_[sequenceOut_];
1473        updateDualsInDual();
1474        abcDj_[sequenceOut_]=0.0;
1475        sequenceOut_=-1;
1476      }
1477      // clear all arrays
1478      clearArrays(-1);
1479    }
1480    // at this stage sequenceIn_ may be <0
1481    if (sequenceIn_<0&&sequenceOut_>=0) {
1482      usefulArray_[arrayForTableauRow_].compact();
1483      // no incoming column is valid
1484      returnCode=noPivotColumn();
1485    }
1486    if (stateOfIteration_==2) {
1487      usefulArray_[arrayForTableauRow_].compact();
1488      sequenceIn_=-1;
1489#ifdef ABC_LONG_FACTORIZATION
1490      abcFactorization_->clearHiddenArrays();
1491#endif
1492      break;
1493    }
1494    if (problemStatus_!=-1) {
1495#ifndef MOVE_UPDATE_WEIGHTS
1496      abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
1497#else
1498      abcDualRowPivot_->updateWeights2(usefulArray_[5],usefulArray_[arrayForFtran_]);
1499#endif
1500#ifdef ABC_LONG_FACTORIZATION
1501      abcFactorization_->clearHiddenArrays();
1502#endif
1503      swapPrimalStuff();
1504      break;
1505    }
1506    if (stateOfIteration_==1) {
1507      // clear all arrays
1508      clearArrays(-2);
1509#ifdef ABC_LONG_FACTORIZATION
1510      abcFactorization_->clearHiddenArrays();
1511#endif
1512    }
1513    if (stateOfIteration_==0) {
1514      // can do these in parallel
1515      // No idea why I need this - but otherwise runs not repeatable (try again??)
1516      //usefulArray_[3].compact();
1517      cilk_spawn updateDualsInDual();
1518      int lastSequenceOut;
1519      int lastDirectionOut;
1520      if (firstFree_<0) {
1521        // can do in parallel
1522        cilk_spawn replaceColumnPart3();
1523        updatePrimalSolution();
1524        swapPrimalStuff();
1525        // dualRow will go to virtual row pivot choice algorithm
1526        // make sure values pass off if it should be
1527        // use Btran array and clear inside dualPivotRow (if used)
1528        lastSequenceOut=sequenceOut_;
1529        lastDirectionOut=directionOut_;
1530        dualPivotRow();
1531        cilk_sync;
1532      } else {
1533        // be more careful as dualPivotRow may do update
1534        cilk_spawn replaceColumnPart3();
1535        updatePrimalSolution();
1536        swapPrimalStuff();
1537        // dualRow will go to virtual row pivot choice algorithm
1538        // make sure values pass off if it should be
1539        // use Btran array and clear inside dualPivotRow (if used)
1540        lastSequenceOut=sequenceOut_;
1541        lastDirectionOut=directionOut_;
1542        cilk_sync;
1543        dualPivotRow();
1544      }
1545      lastPivotRow_=pivotRow_;
1546      if (pivotRow_ >= 0) {
1547        // we found a pivot row
1548        // create as packed
1549        // dual->checkReplacePart1a();
1550        createDualPricingVectorCilk();
1551        swapDualStuff(lastSequenceOut,lastDirectionOut);
1552      }
1553      cilk_sync;
1554    } else {
1555      // after moving dual in values pass
1556      dualPivotRow();
1557      lastPivotRow_=pivotRow_;
1558      if (pivotRow_ >= 0) {
1559        // we found a pivot row
1560        // create as packed
1561        createDualPricingVectorCilk();
1562      }
1563    }
1564    if (pivotRow_<0) {
1565      // no pivot row
1566      clearArrays(-1);
1567      returnCode=noPivotRow();
1568      break;
1569    }
1570    if (numberIterations == numberIterations_&&problemStatus_==-1) {
1571      returnCode=-99;
1572      break;
1573    }
1574  } while (problemStatus_ == -1);
1575  return returnCode;
1576}
1577#if 0
1578void
1579AbcSimplexDual::whileIterating2()
1580{
1581  // get sequenceIn_
1582  dualPivotColumn();
1583  //stateOfIteration_=0;
1584  if (sequenceIn_ >= 0) {
1585    // normal iteration
1586    // update the incoming column
1587    //arrayForReplaceColumn_=getAvailableArray();
1588    getTableauColumnFlipAndStartReplaceCilk();
1589    //usleep(1000);
1590  } else if (stateOfIteration_!=-1) {
1591    stateOfIteration_=2;
1592  }
1593}
1594#endif
1595#endif
1596void 
1597AbcSimplexDual::updatePrimalSolution()
1598{
1599  // finish doing weights
1600#ifndef MOVE_UPDATE_WEIGHTS
1601  abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
1602                                         movement_);
1603#else
1604  abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[5],usefulArray_[arrayForFtran_],
1605                                         movement_);
1606#endif
1607}
1608int xxInfo[6][8];
1609double parallelDual4(AbcSimplexDual * dual)
1610{
1611  int maximumRows=dual->maximumAbcNumberRows();
1612  //int numberTotal=dual->numberTotal();
1613  CoinIndexedVector & update = *dual->usefulArray(dual->arrayForBtran());
1614  CoinPartitionedVector & tableauRow= *dual->usefulArray(dual->arrayForTableauRow());
1615  CoinPartitionedVector & candidateList= *dual->usefulArray(dual->arrayForDualColumn());
1616  AbcMatrix * matrix = dual->abcMatrix();
1617  // do slacks first (move before sync)
1618  double upperTheta=dual->upperTheta();
1619  //assert (upperTheta>-100*dual->dualTolerance()||dual->sequenceIn()>=0||dual->lastFirstFree()>=0);
1620#if ABC_PARALLEL
1621#define NUMBER_BLOCKS NUMBER_ROW_BLOCKS
1622  int numberBlocks=CoinMin(NUMBER_BLOCKS,dual->numberCpus());
1623#else
1624#define NUMBER_BLOCKS 1
1625  int numberBlocks=NUMBER_BLOCKS;
1626#endif
1627#if ABC_PARALLEL
1628  if (dual->parallelMode()==0)
1629    numberBlocks=1;
1630#endif
1631  // see if worth using row copy
1632  bool gotRowCopy = matrix->gotRowCopy();
1633  int number=update.getNumElements();
1634  assert (number);
1635  bool useRowCopy = (gotRowCopy&&(number<<2)<maximumRows);
1636  assert (numberBlocks==matrix->numberRowBlocks());
1637#if ABC_PARALLEL==1
1638  // redo so can pass info in with void *
1639  CoinAbcThreadInfo * infoP = dual->threadInfoPointer();
1640  int cpuMask=((1<<dual->parallelMode())-1);
1641  cpuMask += cpuMask<<5;
1642  dual->setStopStart(cpuMask);
1643#endif
1644  CoinAbcThreadInfo info[NUMBER_BLOCKS];
1645  const int * starts;
1646  if (useRowCopy)
1647    starts=matrix->blockStart();
1648  else
1649    starts=matrix->startColumnBlock();
1650  tableauRow.setPartitions(numberBlocks,starts);
1651  candidateList.setPartitions(numberBlocks,starts);
1652  //printf("free sequence in %d\n",dual->freeSequenceIn());
1653  if (useRowCopy) {
1654    // using row copy
1655#if ABC_PARALLEL
1656    if (numberBlocks>1) {
1657#if ABC_PARALLEL==2
1658    for (int i=0;i<numberBlocks;i++) {
1659      info[i].stuff[1]=i;
1660      info[i].stuff[2]=-1;
1661      info[i].result=upperTheta;
1662      info[i].result= cilk_spawn
1663        matrix->dualColumn1Row(info[i].stuff[1],COIN_DBL_MAX,info[i].stuff[2],
1664                               update,tableauRow,candidateList);
1665    }
1666    cilk_sync;
1667#else
1668      // parallel 1
1669    for (int i=0;i<numberBlocks;i++) {
1670      info[i].status=5;
1671      info[i].stuff[1]=i;
1672      info[i].result=upperTheta;
1673      if (i<numberBlocks-1) {
1674        infoP[i]=info[i];
1675        if (i==numberBlocks-2) 
1676          dual->startParallelStuff(5);
1677      }
1678    }
1679    info[numberBlocks-1].stuff[1]=numberBlocks-1;
1680    info[numberBlocks-1].stuff[2]=-1;
1681    //double freeAlpha;
1682    info[numberBlocks-1].result =
1683      matrix->dualColumn1Row(info[numberBlocks-1].stuff[1],
1684                             upperTheta,info[numberBlocks-1].stuff[2],
1685                             update,tableauRow,candidateList);
1686    dual->stopParallelStuff(5);
1687    for (int i=0;i<numberBlocks-1;i++) 
1688      info[i]=infoP[i];
1689#endif
1690    } else {
1691#endif
1692      info[0].status=5;
1693      info[0].stuff[1]=0;
1694      info[0].result=upperTheta;
1695      info[0].stuff[1]=0;
1696      info[0].stuff[2]=-1;
1697      // worth using row copy
1698      //assert (number>2);
1699      info[0].result =
1700      matrix->dualColumn1Row(info[0].stuff[1],
1701                             upperTheta,info[0].stuff[2],
1702                             update,tableauRow,candidateList);
1703#if ABC_PARALLEL
1704    }
1705#endif
1706  } else { // end of useRowCopy
1707#if ABC_PARALLEL
1708    if (numberBlocks>1) {
1709#if ABC_PARALLEL==2
1710      // do by column
1711      for (int i=0;i<numberBlocks;i++) {
1712        info[i].stuff[1]=i;
1713        info[i].result=upperTheta;
1714        cilk_spawn
1715          matrix->dualColumn1Part(info[i].stuff[1],info[i].stuff[2],
1716                                             info[i].result,
1717                                             update,tableauRow,candidateList);
1718      }
1719      cilk_sync;
1720#else
1721      // parallel 1
1722      // do by column
1723      for (int i=0;i<numberBlocks;i++) {
1724        info[i].status=4;
1725        info[i].stuff[1]=i;
1726        info[i].result=upperTheta;
1727        if (i<numberBlocks-1) {
1728          infoP[i]=info[i];
1729          if (i==numberBlocks-2) 
1730            dual->startParallelStuff(4);
1731        }
1732      }
1733      matrix->dualColumn1Part(info[numberBlocks-1].stuff[1],info[numberBlocks-1].stuff[2],
1734                                                            info[numberBlocks-1].result,
1735                                                            update,tableauRow,candidateList);
1736      dual->stopParallelStuff(4);
1737      for (int i=0;i<numberBlocks-1;i++) 
1738        info[i]=infoP[i];
1739#endif
1740    } else {
1741#endif
1742      info[0].status=4;
1743      info[0].stuff[1]=0;
1744      info[0].result=upperTheta;
1745      info[0].stuff[1]=0;
1746      matrix->dualColumn1Part(info[0].stuff[1],info[0].stuff[2],
1747                                               info[0].result,
1748                                               update,tableauRow,candidateList);
1749#if ABC_PARALLEL
1750    }
1751#endif
1752  }
1753  int sequenceIn[NUMBER_BLOCKS];
1754  bool anyFree=false;
1755#if 0
1756  if (useRowCopy) {
1757    printf("Result at iteration %d slack seqIn %d upperTheta %g\n",
1758           dual->numberIterations(),dual->freeSequenceIn(),upperTheta);
1759    double * arrayT = tableauRow.denseVector();
1760    int * indexT = tableauRow.getIndices();
1761    double * arrayC = candidateList.denseVector();
1762    int * indexC = candidateList.getIndices();
1763    for (int i=0;i<numberBlocks;i++) {
1764      printf("Block %d seqIn %d upperTheta %g first %d last %d firstIn %d nnz %d numrem %d\n",
1765             i,info[i].stuff[2],info[i].result,
1766             xxInfo[0][i],xxInfo[1][i],xxInfo[2][i],xxInfo[3][i],xxInfo[4][i]);
1767      if (xxInfo[3][i]<-35) {
1768        assert (xxInfo[3][i]==tableauRow.getNumElements(i));
1769        assert (xxInfo[4][i]==candidateList.getNumElements(i));
1770        for (int k=0;k<xxInfo[3][i];k++)
1771          printf("T %d %d %g\n",k,indexT[k+xxInfo[2][i]],arrayT[k+xxInfo[2][i]]);
1772        for (int k=0;k<xxInfo[4][i];k++)
1773          printf("C %d %d %g\n",k,indexC[k+xxInfo[2][i]],arrayC[k+xxInfo[2][i]]);
1774      }
1775    }
1776  }
1777#endif
1778  for (int i=0;i<numberBlocks;i++) {
1779    sequenceIn[i]=info[i].stuff[2];
1780    if (sequenceIn[i]>=0)
1781      anyFree=true;
1782    upperTheta=CoinMin(info[i].result,upperTheta);
1783    //assert (info[i].result>-100*dual->dualTolerance()||sequenceIn[i]>=0||dual->lastFirstFree()>=0);
1784  }
1785  if (anyFree) {
1786    int *  COIN_RESTRICT index = tableauRow.getIndices();
1787    double *  COIN_RESTRICT array = tableauRow.denseVector();
1788    // search for free coming in
1789    double bestFree=0.0;
1790    int bestSequence=dual->sequenceIn();
1791    if (bestSequence>=0)
1792      bestFree=dual->alpha();
1793    for (int i=0;i<numberBlocks;i++) {
1794      int iLook=sequenceIn[i];
1795      if (iLook>=0) {
1796        // free variable - search
1797        int start=starts[i];
1798        int end=start+tableauRow.getNumElements(i);
1799        for (int k=start;k<end;k++) {
1800          if (iLook==index[k]) {
1801            if (fabs(bestFree)<fabs(array[k])) {
1802              bestFree=array[k];
1803              bestSequence=iLook;
1804            }
1805            break;
1806          }
1807        }
1808      }
1809    }
1810    if (bestSequence>=0) {
1811      double oldValue = dual->djRegion()[bestSequence];
1812      dual->setSequenceIn(bestSequence);
1813      dual->setAlpha(bestFree);
1814      dual->setTheta(oldValue / bestFree);
1815    }
1816  }
1817  //tableauRow.compact();
1818  //candidateList.compact();
1819#if 0//ndef NDEBUG
1820  tableauRow.setPackedMode(true);
1821  tableauRow.sortPacked();
1822  candidateList.setPackedMode(true);
1823  candidateList.sortPacked();
1824#endif
1825  return upperTheta;
1826}
1827#if ABC_PARALLEL==2
1828static void
1829parallelDual5a(AbcSimplexFactorization * factorization,
1830              CoinIndexedVector * whichVector,
1831              int numberCpu,
1832              int whichCpu,
1833              double * weights)
1834{
1835  double *  COIN_RESTRICT array = whichVector->denseVector();
1836  int *  COIN_RESTRICT which = whichVector->getIndices();
1837  int numberRows=factorization->numberRows();
1838  for (int i = whichCpu; i < numberRows; i+=numberCpu) {
1839    double value = 0.0;
1840    array[i] = 1.0;
1841    which[0] = i;
1842    whichVector->setNumElements(1);
1843    factorization->updateColumnTransposeCpu(*whichVector,whichCpu);
1844    int number = whichVector->getNumElements();
1845    for (int j = 0; j < number; j++) {
1846      int k= which[j];
1847      value += array[k] * array[k];
1848      array[k] = 0.0;
1849    }
1850    weights[i] = value;
1851  }
1852  whichVector->setNumElements(0);
1853}
1854#endif
1855#if ABC_PARALLEL==2
1856void
1857parallelDual5(AbcSimplexFactorization * factorization,
1858              CoinIndexedVector ** whichVector,
1859              int numberCpu,
1860              int whichCpu,
1861              double * weights)
1862{
1863  if (whichCpu) {
1864    cilk_spawn parallelDual5(factorization,whichVector,numberCpu,whichCpu-1,weights);
1865    parallelDual5a(factorization,whichVector[whichCpu],numberCpu,whichCpu,weights);
1866    cilk_sync;
1867  } else {
1868    parallelDual5a(factorization,whichVector[whichCpu],numberCpu,whichCpu,weights);
1869  }
1870}
1871#endif
1872// cilk seems a bit fragile
1873#define CILK_FRAGILE 1
1874#if CILK_FRAGILE>1
1875#undef cilk_spawn
1876#undef cilk_sync
1877#define cilk_spawn
1878#define cilk_sync
1879#define ONWARD 0
1880#elif CILK_FRAGILE==1
1881#define ONWARD 0
1882#else
1883#define ONWARD 1
1884#endif
1885// Code below is just a translation of LAPACK
1886#define BLOCKING1 8 // factorization strip
1887#define BLOCKING2 8 // dgemm recursive
1888#define BLOCKING3 16 // dgemm parallel
1889/* type
1890   0 Left Lower NoTranspose Unit
1891   1 Left Upper NoTranspose NonUnit
1892   2 Left Lower Transpose Unit
1893   3 Left Upper Transpose NonUnit
1894*/
1895static void CoinAbcDtrsmFactor(int m, int n, double * COIN_RESTRICT a,int lda)
1896{
1897  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
1898  assert (m==BLOCKING8);
1899  // 0 Left Lower NoTranspose Unit
1900  /* entry for column j and row i (when multiple of BLOCKING8)
1901     is at aBlocked+j*m+i*BLOCKING8
1902  */
1903  double * COIN_RESTRICT aBase2 = a;
1904  double * COIN_RESTRICT bBase2 = aBase2+lda*BLOCKING8;
1905  for (int jj=0;jj<n;jj+=BLOCKING8) {
1906    double * COIN_RESTRICT bBase = bBase2;
1907    for (int j=jj;j<jj+BLOCKING8;j++) {
1908#if 0
1909      double * COIN_RESTRICT aBase = aBase2;
1910      for (int k=0;k<BLOCKING8;k++) {
1911        double bValue = bBase[k];
1912        if (bValue) {
1913          for (int i=k+1;i<BLOCKING8;i++) {
1914            bBase[i]-=bValue*aBase[i];
1915          }
1916        }
1917        aBase+=BLOCKING8;
1918      }
1919#else
1920      // unroll - stay in registers - don't test for zeros
1921      assert (BLOCKING8==8);
1922      double bValue0 = bBase[0];
1923      double bValue1 = bBase[1];
1924      double bValue2 = bBase[2];
1925      double bValue3 = bBase[3];
1926      double bValue4 = bBase[4];
1927      double bValue5 = bBase[5];
1928      double bValue6 = bBase[6];
1929      double bValue7 = bBase[7];
1930      bValue1-=bValue0*a[1+0*BLOCKING8];
1931      bBase[1] = bValue1;
1932      bValue2-=bValue0*a[2+0*BLOCKING8];
1933      bValue3-=bValue0*a[3+0*BLOCKING8];
1934      bValue4-=bValue0*a[4+0*BLOCKING8];
1935      bValue5-=bValue0*a[5+0*BLOCKING8];
1936      bValue6-=bValue0*a[6+0*BLOCKING8];
1937      bValue7-=bValue0*a[7+0*BLOCKING8];
1938      bValue2-=bValue1*a[2+1*BLOCKING8];
1939      bBase[2] = bValue2;
1940      bValue3-=bValue1*a[3+1*BLOCKING8];
1941      bValue4-=bValue1*a[4+1*BLOCKING8];
1942      bValue5-=bValue1*a[5+1*BLOCKING8];
1943      bValue6-=bValue1*a[6+1*BLOCKING8];
1944      bValue7-=bValue1*a[7+1*BLOCKING8];
1945      bValue3-=bValue2*a[3+2*BLOCKING8];
1946      bBase[3] = bValue3;
1947      bValue4-=bValue2*a[4+2*BLOCKING8];
1948      bValue5-=bValue2*a[5+2*BLOCKING8];
1949      bValue6-=bValue2*a[6+2*BLOCKING8];
1950      bValue7-=bValue2*a[7+2*BLOCKING8];
1951      bValue4-=bValue3*a[4+3*BLOCKING8];
1952      bBase[4] = bValue4;
1953      bValue5-=bValue3*a[5+3*BLOCKING8];
1954      bValue6-=bValue3*a[6+3*BLOCKING8];
1955      bValue7-=bValue3*a[7+3*BLOCKING8];
1956      bValue5-=bValue4*a[5+4*BLOCKING8];
1957      bBase[5] = bValue5;
1958      bValue6-=bValue4*a[6+4*BLOCKING8];
1959      bValue7-=bValue4*a[7+4*BLOCKING8];
1960      bValue6-=bValue5*a[6+5*BLOCKING8];
1961      bBase[6] = bValue6;
1962      bValue7-=bValue5*a[7+5*BLOCKING8];
1963      bValue7-=bValue6*a[7+6*BLOCKING8];
1964      bBase[7] = bValue7;
1965#endif
1966      bBase += BLOCKING8;
1967    }
1968    bBase2 +=lda*BLOCKING8; 
1969  }
1970}
1971#define UNROLL_DTRSM 16
1972#define CILK_DTRSM 32
1973static void dtrsm0(int kkk, int first, int last,
1974                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
1975{
1976  int mm=CoinMin(kkk+UNROLL_DTRSM*BLOCKING8,m);
1977  assert ((last-first)%BLOCKING8==0);
1978  if (last-first>CILK_DTRSM) {
1979    int mid=((first+last)>>4)<<3;
1980    cilk_spawn dtrsm0(kkk,first,mid,m,a,b);
1981    dtrsm0(kkk,mid,last,m,a,b);
1982    cilk_sync;
1983  } else {
1984    const double * COIN_RESTRICT aBaseA = a+UNROLL_DTRSM*BLOCKING8X8+kkk*BLOCKING8;
1985    aBaseA += (first-mm)*BLOCKING8-BLOCKING8X8;
1986    aBaseA += m*kkk;
1987    for (int ii=first;ii<last;ii+=BLOCKING8) {
1988      aBaseA += BLOCKING8X8;
1989      const double * COIN_RESTRICT aBaseB=aBaseA;
1990      double * COIN_RESTRICT bBaseI = b+ii;
1991      for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
1992        double * COIN_RESTRICT bBase = b+kk;
1993        const double * COIN_RESTRICT aBase2 = aBaseB;//a+UNROLL_DTRSM*BLOCKING8X8+m*kk+kkk*BLOCKING8;
1994        //aBase2 += (ii-mm)*BLOCKING8;
1995        //assert (aBase2==aBaseB);
1996        aBaseB+=m*BLOCKING8;
1997#if AVX2 !=2
1998#define ALTERNATE_INNER
1999#ifndef ALTERNATE_INNER
2000        for (int k=0;k<BLOCKING8;k++) {
2001          //coin_prefetch_const(aBase2+BLOCKING8);             
2002          for (int i=0;i<BLOCKING8;i++) {
2003            bBaseI[i] -= bBase[k]*aBase2[i];
2004          }
2005          aBase2 += BLOCKING8;
2006        }
2007#else
2008        double b0=bBaseI[0];
2009        double b1=bBaseI[1];
2010        double b2=bBaseI[2];
2011        double b3=bBaseI[3];
2012        double b4=bBaseI[4];
2013        double b5=bBaseI[5];
2014        double b6=bBaseI[6];
2015        double b7=bBaseI[7];
2016        for (int k=0;k<BLOCKING8;k++) {
2017          //coin_prefetch_const(aBase2+BLOCKING8);             
2018          double bValue=bBase[k];
2019          b0 -= bValue * aBase2[0];
2020          b1 -= bValue * aBase2[1];
2021          b2 -= bValue * aBase2[2];
2022          b3 -= bValue * aBase2[3];
2023          b4 -= bValue * aBase2[4];
2024          b5 -= bValue * aBase2[5];
2025          b6 -= bValue * aBase2[6];
2026          b7 -= bValue * aBase2[7];
2027          aBase2 += BLOCKING8;
2028        }
2029        bBaseI[0]=b0;
2030        bBaseI[1]=b1;
2031        bBaseI[2]=b2;
2032        bBaseI[3]=b3;
2033        bBaseI[4]=b4;
2034        bBaseI[5]=b5;
2035        bBaseI[6]=b6;
2036        bBaseI[7]=b7;
2037#endif
2038#else
2039        __m256d b0=_mm256_load_pd(bBaseI);
2040        __m256d b1=_mm256_load_pd(bBaseI+4);
2041        for (int k=0;k<BLOCKING8;k++) {
2042          //coin_prefetch_const(aBase2+BLOCKING8);             
2043          __m256d bb = _mm256_broadcast_sd(bBase+k);
2044          __m256d a0 = _mm256_load_pd(aBase2);
2045          b0 -= a0*bb;
2046          __m256d a1 = _mm256_load_pd(aBase2+4);
2047          b1 -= a1*bb;
2048          aBase2 += BLOCKING8;
2049        }
2050        //_mm256_store_pd ((bBaseI), (b0));
2051        *reinterpret_cast<__m256d *>(bBaseI)=b0;
2052        //_mm256_store_pd (bBaseI+4, b1);
2053        *reinterpret_cast<__m256d *>(bBaseI+4)=b1;
2054#endif
2055      }
2056    }
2057  }
2058}
2059void CoinAbcDtrsm0(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2060{
2061  assert ((m&(BLOCKING8-1))==0);
2062  // 0 Left Lower NoTranspose Unit
2063  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM*BLOCKING8) {
2064    int mm=CoinMin(m,kkk+UNROLL_DTRSM*BLOCKING8);
2065    for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
2066      const double * COIN_RESTRICT aBase2 = a+kk*(m+BLOCKING8);
2067      double * COIN_RESTRICT bBase = b+kk;
2068      for (int k=0;k<BLOCKING8;k++) {
2069        for (int i=k+1;i<BLOCKING8;i++) {
2070          bBase[i] -= bBase[k]*aBase2[i];
2071        }
2072        aBase2 += BLOCKING8;
2073      }
2074      for (int ii=kk+BLOCKING8;ii<mm;ii+=BLOCKING8) {
2075        double * COIN_RESTRICT bBaseI = b+ii;
2076        for (int k=0;k<BLOCKING8;k++) {
2077          //coin_prefetch_const(aBase2+BLOCKING8);             
2078          for (int i=0;i<BLOCKING8;i++) {
2079            bBaseI[i] -= bBase[k]*aBase2[i];
2080          }
2081          aBase2 += BLOCKING8;
2082        }
2083      }
2084    }
2085    dtrsm0(kkk,mm,m,m,a,b);
2086  }
2087}
2088static void dtrsm1(int kkk, int first, int last,
2089                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2090{
2091  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2092  assert ((last-first)%BLOCKING8==0);
2093  if (last-first>CILK_DTRSM) {
2094    int mid=((first+last)>>4)<<3;
2095    cilk_spawn dtrsm1(kkk,first,mid,m,a,b);
2096    dtrsm1(kkk,mid,last,m,a,b);
2097    cilk_sync;
2098  } else {
2099    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2100      double * COIN_RESTRICT bBase2 = b+iii;
2101      const double * COIN_RESTRICT aBaseA=
2102        a+BLOCKING8X8+BLOCKING8*iii;
2103      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2104        double * COIN_RESTRICT bBase = b+ii;
2105        const double * COIN_RESTRICT aBase=aBaseA+m*ii;
2106#if AVX2 !=2
2107#ifndef ALTERNATE_INNER
2108        for (int i=BLOCKING8-1;i>=0;i--) {
2109          aBase -= BLOCKING8;
2110          //coin_prefetch_const(aBase-BLOCKING8);               
2111          for (int k=BLOCKING8-1;k>=0;k--) {
2112            bBase2[k] -= bBase[i]*aBase[k];
2113          }
2114        }
2115#else
2116        double b0=bBase2[0];
2117        double b1=bBase2[1];
2118        double b2=bBase2[2];
2119        double b3=bBase2[3];
2120        double b4=bBase2[4];
2121        double b5=bBase2[5];
2122        double b6=bBase2[6];
2123        double b7=bBase2[7];
2124        for (int i=BLOCKING8-1;i>=0;i--) {
2125          aBase -= BLOCKING8;
2126          //coin_prefetch_const(aBase-BLOCKING8);               
2127          double bValue=bBase[i];
2128          b0 -= bValue * aBase[0];
2129          b1 -= bValue * aBase[1];
2130          b2 -= bValue * aBase[2];
2131          b3 -= bValue * aBase[3];
2132          b4 -= bValue * aBase[4];
2133          b5 -= bValue * aBase[5];
2134          b6 -= bValue * aBase[6];
2135          b7 -= bValue * aBase[7];
2136        }
2137        bBase2[0]=b0;
2138        bBase2[1]=b1;
2139        bBase2[2]=b2;
2140        bBase2[3]=b3;
2141        bBase2[4]=b4;
2142        bBase2[5]=b5;
2143        bBase2[6]=b6;
2144        bBase2[7]=b7;
2145#endif
2146#else
2147        __m256d b0=_mm256_load_pd(bBase2);
2148        __m256d b1=_mm256_load_pd(bBase2+4);
2149        for (int i=BLOCKING8-1;i>=0;i--) {
2150          aBase -= BLOCKING8;
2151          //coin_prefetch_const(aBase5-BLOCKING8);             
2152          __m256d bb = _mm256_broadcast_sd(bBase+i);
2153          __m256d a0 = _mm256_load_pd(aBase);
2154          b0 -= a0*bb;
2155          __m256d a1 = _mm256_load_pd(aBase+4);
2156          b1 -= a1*bb;
2157        }
2158        //_mm256_store_pd (bBase2, b0);
2159        *reinterpret_cast<__m256d *>(bBase2)=b0;
2160        //_mm256_store_pd (bBase2+4, b1);
2161        *reinterpret_cast<__m256d *>(bBase2+4)=b1;
2162#endif
2163      }
2164    }
2165  }
2166}
2167void CoinAbcDtrsm1(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2168{
2169  assert ((m&(BLOCKING8-1))==0);
2170  // 1 Left Upper NoTranspose NonUnit
2171  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2172    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2173    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2174      const double * COIN_RESTRICT aBase = a+m*ii+ii*BLOCKING8;
2175      double * COIN_RESTRICT bBase = b+ii;
2176      for (int i=BLOCKING8-1;i>=0;i--) {
2177        bBase[i] *= aBase[i*(BLOCKING8+1)];
2178        for (int k=i-1;k>=0;k--) {
2179          bBase[k] -= bBase[i]*aBase[k+i*BLOCKING8];
2180        }
2181      }
2182      double * COIN_RESTRICT bBase2 = bBase;
2183      for (int iii=ii;iii>mm;iii-=BLOCKING8) {
2184        bBase2 -= BLOCKING8;
2185        for (int i=BLOCKING8-1;i>=0;i--) {
2186          aBase -= BLOCKING8;
2187          coin_prefetch_const(aBase-BLOCKING8);         
2188          for (int k=BLOCKING8-1;k>=0;k--) {
2189            bBase2[k] -= bBase[i]*aBase[k];
2190          }
2191        }
2192      }
2193    }
2194    dtrsm1(kkk, 0, mm,  m, a, b);
2195  }
2196}
2197static void dtrsm2(int kkk, int first, int last,
2198                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2199{
2200  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2201  assert ((last-first)%BLOCKING8==0);
2202  if (last-first>CILK_DTRSM) {
2203    int mid=((first+last)>>4)<<3;
2204    cilk_spawn dtrsm2(kkk,first,mid,m,a,b);
2205    dtrsm2(kkk,mid,last,m,a,b);
2206    cilk_sync;
2207  } else {
2208    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2209      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2210        const double * COIN_RESTRICT bBase = b+ii;
2211        double * COIN_RESTRICT bBase2=b+iii;
2212        const double * COIN_RESTRICT aBase=a+ii*BLOCKING8+m*iii+BLOCKING8X8;
2213        for (int j=BLOCKING8-1;j>=0;j-=4) {
2214          double bValue0=bBase2[j];
2215          double bValue1=bBase2[j-1];
2216          double bValue2=bBase2[j-2];
2217          double bValue3=bBase2[j-3];
2218          bValue0 -= aBase[-1*BLOCKING8+7]*bBase[7];
2219          bValue1 -= aBase[-2*BLOCKING8+7]*bBase[7];
2220          bValue2 -= aBase[-3*BLOCKING8+7]*bBase[7];
2221          bValue3 -= aBase[-4*BLOCKING8+7]*bBase[7];
2222          bValue0 -= aBase[-1*BLOCKING8+6]*bBase[6];
2223          bValue1 -= aBase[-2*BLOCKING8+6]*bBase[6];
2224          bValue2 -= aBase[-3*BLOCKING8+6]*bBase[6];
2225          bValue3 -= aBase[-4*BLOCKING8+6]*bBase[6];
2226          bValue0 -= aBase[-1*BLOCKING8+5]*bBase[5];
2227          bValue1 -= aBase[-2*BLOCKING8+5]*bBase[5];
2228          bValue2 -= aBase[-3*BLOCKING8+5]*bBase[5];
2229          bValue3 -= aBase[-4*BLOCKING8+5]*bBase[5];
2230          bValue0 -= aBase[-1*BLOCKING8+4]*bBase[4];
2231          bValue1 -= aBase[-2*BLOCKING8+4]*bBase[4];
2232          bValue2 -= aBase[-3*BLOCKING8+4]*bBase[4];
2233          bValue3 -= aBase[-4*BLOCKING8+4]*bBase[4];
2234          bValue0 -= aBase[-1*BLOCKING8+3]*bBase[3];
2235          bValue1 -= aBase[-2*BLOCKING8+3]*bBase[3];
2236          bValue2 -= aBase[-3*BLOCKING8+3]*bBase[3];
2237          bValue3 -= aBase[-4*BLOCKING8+3]*bBase[3];
2238          bValue0 -= aBase[-1*BLOCKING8+2]*bBase[2];
2239          bValue1 -= aBase[-2*BLOCKING8+2]*bBase[2];
2240          bValue2 -= aBase[-3*BLOCKING8+2]*bBase[2];
2241          bValue3 -= aBase[-4*BLOCKING8+2]*bBase[2];
2242          bValue0 -= aBase[-1*BLOCKING8+1]*bBase[1];
2243          bValue1 -= aBase[-2*BLOCKING8+1]*bBase[1];
2244          bValue2 -= aBase[-3*BLOCKING8+1]*bBase[1];
2245          bValue3 -= aBase[-4*BLOCKING8+1]*bBase[1];
2246          bValue0 -= aBase[-1*BLOCKING8+0]*bBase[0];
2247          bValue1 -= aBase[-2*BLOCKING8+0]*bBase[0];
2248          bValue2 -= aBase[-3*BLOCKING8+0]*bBase[0];
2249          bValue3 -= aBase[-4*BLOCKING8+0]*bBase[0];
2250          bBase2[j]=bValue0;
2251          bBase2[j-1]=bValue1;
2252          bBase2[j-2]=bValue2;
2253          bBase2[j-3]=bValue3;
2254          aBase -= 4*BLOCKING8;
2255        }
2256      }
2257    }
2258  }
2259}
2260void CoinAbcDtrsm2(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2261{
2262  assert ((m&(BLOCKING8-1))==0);
2263  // 2 Left Lower Transpose Unit
2264  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2265    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2266    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2267      double * COIN_RESTRICT bBase = b+ii;
2268      double * COIN_RESTRICT bBase2 = bBase;
2269      const double * COIN_RESTRICT aBase=a+m*ii+(ii+BLOCKING8)*BLOCKING8;
2270      for (int i=BLOCKING8-1;i>=0;i--) {
2271        aBase -= BLOCKING8;
2272        for (int k=i+1;k<BLOCKING8;k++) {
2273          bBase2[i] -= aBase[k]*bBase2[k];
2274        }
2275      }
2276      for (int iii=ii-BLOCKING8;iii>=mm;iii-=BLOCKING8) {
2277        bBase2 -= BLOCKING8;
2278        assert (bBase2==b+iii);
2279        aBase -= m*BLOCKING8;
2280        const double * COIN_RESTRICT aBase2 = aBase+BLOCKING8X8;
2281        for (int i=BLOCKING8-1;i>=0;i--) {
2282          aBase2 -= BLOCKING8;
2283          for (int k=BLOCKING8-1;k>=0;k--) {
2284            bBase2[i] -= aBase2[k]*bBase[k];
2285          }
2286        }
2287      }
2288    }
2289    dtrsm2(kkk, 0, mm,  m, a, b);
2290  }
2291}
2292#define UNROLL_DTRSM3 16
2293#define CILK_DTRSM3 32
2294static void dtrsm3(int kkk, int first, int last,
2295                   int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2296{
2297  //int mm=CoinMin(kkk+UNROLL_DTRSM3*BLOCKING8,m);
2298  assert ((last-first)%BLOCKING8==0);
2299  if (last-first>CILK_DTRSM3) {
2300    int mid=((first+last)>>4)<<3;
2301    cilk_spawn dtrsm3(kkk,first,mid,m,a,b);
2302    dtrsm3(kkk,mid,last,m,a,b);
2303    cilk_sync;
2304  } else {
2305    for (int kk=0;kk<kkk;kk+=BLOCKING8) {
2306      for (int ii=first;ii<last;ii+=BLOCKING8) {
2307        double * COIN_RESTRICT bBase = b+ii;
2308        double * COIN_RESTRICT bBase2 = b+kk;
2309        const double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
2310        for (int j=0;j<BLOCKING8;j+=4) {
2311          double bValue0=bBase[j];
2312          double bValue1=bBase[j+1];
2313          double bValue2=bBase[j+2];
2314          double bValue3=bBase[j+3];
2315          bValue0 -= aBase[0]*bBase2[0];
2316          bValue1 -= aBase[1*BLOCKING8+0]*bBase2[0];
2317          bValue2 -= aBase[2*BLOCKING8+0]*bBase2[0];
2318          bValue3 -= aBase[3*BLOCKING8+0]*bBase2[0];
2319          bValue0 -= aBase[1]*bBase2[1];
2320          bValue1 -= aBase[1*BLOCKING8+1]*bBase2[1];
2321          bValue2 -= aBase[2*BLOCKING8+1]*bBase2[1];
2322          bValue3 -= aBase[3*BLOCKING8+1]*bBase2[1];
2323          bValue0 -= aBase[2]*bBase2[2];
2324          bValue1 -= aBase[1*BLOCKING8+2]*bBase2[2];
2325          bValue2 -= aBase[2*BLOCKING8+2]*bBase2[2];
2326          bValue3 -= aBase[3*BLOCKING8+2]*bBase2[2];
2327          bValue0 -= aBase[3]*bBase2[3];
2328          bValue1 -= aBase[1*BLOCKING8+3]*bBase2[3];
2329          bValue2 -= aBase[2*BLOCKING8+3]*bBase2[3];
2330          bValue3 -= aBase[3*BLOCKING8+3]*bBase2[3];
2331          bValue0 -= aBase[4]*bBase2[4];
2332          bValue1 -= aBase[1*BLOCKING8+4]*bBase2[4];
2333          bValue2 -= aBase[2*BLOCKING8+4]*bBase2[4];
2334          bValue3 -= aBase[3*BLOCKING8+4]*bBase2[4];
2335          bValue0 -= aBase[5]*bBase2[5];
2336          bValue1 -= aBase[1*BLOCKING8+5]*bBase2[5];
2337          bValue2 -= aBase[2*BLOCKING8+5]*bBase2[5];
2338          bValue3 -= aBase[3*BLOCKING8+5]*bBase2[5];
2339          bValue0 -= aBase[6]*bBase2[6];
2340          bValue1 -= aBase[1*BLOCKING8+6]*bBase2[6];
2341          bValue2 -= aBase[2*BLOCKING8+7]*bBase2[7];
2342          bValue3 -= aBase[3*BLOCKING8+6]*bBase2[6];
2343          bValue0 -= aBase[7]*bBase2[7];
2344          bValue1 -= aBase[1*BLOCKING8+7]*bBase2[7];
2345          bValue2 -= aBase[2*BLOCKING8+6]*bBase2[6];
2346          bValue3 -= aBase[3*BLOCKING8+7]*bBase2[7];
2347          bBase[j]=bValue0;
2348          bBase[j+1]=bValue1;
2349          bBase[j+2]=bValue2;
2350          bBase[j+3]=bValue3;
2351          aBase += 4*BLOCKING8;
2352        }
2353      }
2354    }
2355  }
2356}
2357void CoinAbcDtrsm3(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
2358{
2359  assert ((m&(BLOCKING8-1))==0);
2360  // 3 Left Upper Transpose NonUnit
2361  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM3*BLOCKING8) {
2362    int mm=CoinMin(m,kkk+UNROLL_DTRSM3*BLOCKING8);
2363    dtrsm3(kkk, kkk, mm,  m, a, b);
2364    for (int ii=kkk;ii<mm;ii+=BLOCKING8) {
2365      double * COIN_RESTRICT bBase = b+ii;
2366      for (int kk=kkk;kk<ii;kk+=BLOCKING8) {
2367        double * COIN_RESTRICT bBase2 = b+kk;
2368        const double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
2369        for (int i=0;i<BLOCKING8;i++) {
2370          for (int k=0;k<BLOCKING8;k++) {
2371            bBase[i] -= aBase[k]*bBase2[k];
2372          }
2373          aBase += BLOCKING8;
2374        }
2375      }
2376      const double * COIN_RESTRICT aBase=a+ii*m+ii*BLOCKING8;
2377      for (int i=0;i<BLOCKING8;i++) {
2378        for (int k=0;k<i;k++) {
2379          bBase[i] -= aBase[k]*bBase[k];
2380        }
2381        bBase[i] = bBase[i]*aBase[i];
2382        aBase += BLOCKING8;
2383      }
2384    }
2385  }
2386}
2387static void 
2388CoinAbcDlaswp(int n, double * COIN_RESTRICT a,int lda,int start,int end, int * ipiv)
2389{
2390  assert ((n&(BLOCKING8-1))==0);
2391  double * COIN_RESTRICT aThis3 = a;
2392  for (int j=0;j<n;j+=BLOCKING8) {
2393    for (int i=start;i<end;i++) {
2394      int ip=ipiv[i];
2395      if (ip!=i) {
2396        double * COIN_RESTRICT aTop=aThis3+j*lda;
2397        int iBias = ip/BLOCKING8;
2398        ip-=iBias*BLOCKING8;
2399        double * COIN_RESTRICT aNotTop=aTop+iBias*BLOCKING8X8;
2400        iBias = i/BLOCKING8;
2401        int i2=i-iBias*BLOCKING8;
2402        aTop += iBias*BLOCKING8X8;
2403        for (int k=0;k<BLOCKING8;k++) {
2404          double temp=aTop[i2];
2405          aTop[i2]=aNotTop[ip];
2406          aNotTop[ip]=temp;
2407          aTop += BLOCKING8;
2408          aNotTop += BLOCKING8;
2409        }
2410      }
2411    }
2412  }
2413}
2414extern void CoinAbcDgemm(int m, int n, int k, double * COIN_RESTRICT a,int lda,
2415                          double * COIN_RESTRICT b,double * COIN_RESTRICT c
2416#if ABC_PARALLEL==2
2417                          ,int parallelMode
2418#endif
2419                         );
2420static int CoinAbcDgetrf2(int m, int n, double * COIN_RESTRICT a, int * ipiv)
2421{
2422  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
2423  assert (n==BLOCKING8);
2424  int dimension=CoinMin(m,n);
2425  /* entry for column j and row i (when multiple of BLOCKING8)
2426     is at aBlocked+j*m+i*BLOCKING8
2427  */
2428  assert (dimension==BLOCKING8);
2429  //printf("Dgetrf2 m=%d n=%d\n",m,n);
2430  double * COIN_RESTRICT aThis3 = a;
2431  for (int j=0;j<dimension;j++) {
2432    int pivotRow=-1;
2433    double largest=0.0;
2434    double * COIN_RESTRICT aThis2 = aThis3+j*BLOCKING8;
2435    // this block
2436    int pivotRow2=-1;
2437    for (int i=j;i<BLOCKING8;i++) {
2438      double value=fabs(aThis2[i]);
2439      if (value>largest) {
2440        largest=value;
2441        pivotRow2=i;
2442      }
2443    }
2444    // other blocks
2445    int iBias=0;
2446    for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
2447      aThis2+=BLOCKING8*BLOCKING8;
2448      for (int i=0;i<BLOCKING8;i++) {
2449        double value=fabs(aThis2[i]);
2450        if (value>largest) {
2451          largest=value;
2452          pivotRow2=i;
2453          iBias=ii;
2454        }
2455      }
2456    }
2457    pivotRow=pivotRow2+iBias;
2458    ipiv[j]=pivotRow;
2459    if (largest) {
2460      if (j!=pivotRow) {
2461        double * COIN_RESTRICT aTop=aThis3;
2462        double * COIN_RESTRICT aNotTop=aThis3+iBias*BLOCKING8;
2463        for (int i=0;i<n;i++) {
2464          double value = aTop[j];
2465          aTop[j]=aNotTop[pivotRow2];
2466          aNotTop[pivotRow2]=value;
2467          aTop += BLOCKING8;
2468          aNotTop += BLOCKING8;
2469        }
2470      }
2471      aThis2 = aThis3+j*BLOCKING8;
2472      double pivotMultiplier = 1.0 / aThis2[j];
2473      aThis2[j] = pivotMultiplier;
2474      // Do L
2475      // this block
2476      for (int i=j+1;i<BLOCKING8;i++) 
2477        aThis2[i] *= pivotMultiplier;
2478      // other blocks
2479      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
2480        aThis2+=BLOCKING8*BLOCKING8;
2481        for (int i=0;i<BLOCKING8;i++) 
2482          aThis2[i] *= pivotMultiplier;
2483      }
2484      // update rest
2485      double * COIN_RESTRICT aPut;
2486      aThis2 = aThis3+j*BLOCKING8;
2487      aPut = aThis2+BLOCKING8;
2488      double * COIN_RESTRICT aPut2 = aPut;
2489      // this block
2490      for (int i=j+1;i<BLOCKING8;i++) {
2491        double value = aPut[j];
2492        for (int k=j+1;k<BLOCKING8;k++) 
2493          aPut[k] -= value*aThis2[k];
2494        aPut += BLOCKING8;
2495      }
2496      // other blocks
2497      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
2498        aThis2 += BLOCKING8*BLOCKING8;
2499        aPut = aThis2+BLOCKING8;
2500        double * COIN_RESTRICT aPut2a = aPut2;
2501        for (int i=j+1;i<BLOCKING8;i++) {
2502          double value = aPut2a[j];
2503          for (int k=0;k<BLOCKING8;k++) 
2504            aPut[k] -= value*aThis2[k];
2505          aPut += BLOCKING8;
2506          aPut2a += BLOCKING8;
2507        }
2508      }
2509    } else {
2510      return 1;
2511    }
2512  }
2513  return 0;
2514}
2515
2516int 
2517CoinAbcDgetrf(int m, int n, double * COIN_RESTRICT a, int lda, int * ipiv
2518#if ABC_PARALLEL==2
2519                          ,int parallelMode
2520#endif
2521)
2522{
2523  assert (m==n);
2524  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
2525  if (m<BLOCKING8) {
2526    return CoinAbcDgetrf2(m,n,a,ipiv);
2527  } else {
2528    for (int j=0;j<n;j+=BLOCKING8) {
2529      int start=j;
2530      int newSize=CoinMin(BLOCKING8,n-j);
2531      int end=j+newSize;
2532      int returnCode=CoinAbcDgetrf2(m-start,newSize,a+(start*lda+start*BLOCKING8),
2533                                     ipiv+start);
2534      if (!returnCode) {
2535        // adjust
2536        for (int k=start;k<end;k++)
2537          ipiv[k]+=start;
2538        // swap 0<start
2539        CoinAbcDlaswp(start,a,lda,start,end,ipiv);
2540        if (end<n) {
2541          // swap >=end
2542          CoinAbcDlaswp(n-end,a+end*lda,lda,start,end,ipiv);
2543          CoinAbcDtrsmFactor(newSize,n-end,a+(start*lda+start*BLOCKING8),lda);
2544          CoinAbcDgemm(n-end,n-end,newSize,
2545                        a+start*lda+end*BLOCKING8,lda,
2546                        a+end*lda+start*BLOCKING8,a+end*lda+end*BLOCKING8
2547#if ABC_PARALLEL==2
2548                          ,parallelMode
2549#endif
2550                        );
2551        }
2552      } else {
2553        return returnCode;
2554      }
2555    }
2556  }
2557  return 0;
2558}
2559void
2560CoinAbcDgetrs(char trans,int m, double * COIN_RESTRICT a, double * COIN_RESTRICT work)
2561{
2562  assert ((m&(BLOCKING8-1))==0);
2563  if (trans=='N') {
2564    //CoinAbcDlaswp1(work,m,ipiv);
2565    CoinAbcDtrsm0(m,a,work);
2566    CoinAbcDtrsm1(m,a,work);
2567  } else {
2568    assert (trans=='T');
2569    CoinAbcDtrsm3(m,a,work);
2570    CoinAbcDtrsm2(m,a,work);
2571    //CoinAbcDlaswp1Backwards(work,m,ipiv);
2572  }
2573}
2574#ifdef ABC_LONG_FACTORIZATION
2575/// ****** Start long double version
2576/* type
2577   0 Left Lower NoTranspose Unit
2578   1 Left Upper NoTranspose NonUnit
2579   2 Left Lower Transpose Unit
2580   3 Left Upper Transpose NonUnit
2581*/
2582static void CoinAbcDtrsmFactor(int m, int n, long double * COIN_RESTRICT a,int lda)
2583{
2584  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
2585  assert (m==BLOCKING8);
2586  // 0 Left Lower NoTranspose Unit
2587  /* entry for column j and row i (when multiple of BLOCKING8)
2588     is at aBlocked+j*m+i*BLOCKING8
2589  */
2590  long double * COIN_RESTRICT aBase2 = a;
2591  long double * COIN_RESTRICT bBase2 = aBase2+lda*BLOCKING8;
2592  for (int jj=0;jj<n;jj+=BLOCKING8) {
2593    long double * COIN_RESTRICT bBase = bBase2;
2594    for (int j=jj;j<jj+BLOCKING8;j++) {
2595#if 0
2596      long double * COIN_RESTRICT aBase = aBase2;
2597      for (int k=0;k<BLOCKING8;k++) {
2598        long double bValue = bBase[k];
2599        if (bValue) {
2600          for (int i=k+1;i<BLOCKING8;i++) {
2601            bBase[i]-=bValue*aBase[i];
2602          }
2603        }
2604        aBase+=BLOCKING8;
2605      }
2606#else
2607      // unroll - stay in registers - don't test for zeros
2608      assert (BLOCKING8==8);
2609      long double bValue0 = bBase[0];
2610      long double bValue1 = bBase[1];
2611      long double bValue2 = bBase[2];
2612      long double bValue3 = bBase[3];
2613      long double bValue4 = bBase[4];
2614      long double bValue5 = bBase[5];
2615      long double bValue6 = bBase[6];
2616      long double bValue7 = bBase[7];
2617      bValue1-=bValue0*a[1+0*BLOCKING8];
2618      bBase[1] = bValue1;
2619      bValue2-=bValue0*a[2+0*BLOCKING8];
2620      bValue3-=bValue0*a[3+0*BLOCKING8];
2621      bValue4-=bValue0*a[4+0*BLOCKING8];
2622      bValue5-=bValue0*a[5+0*BLOCKING8];
2623      bValue6-=bValue0*a[6+0*BLOCKING8];
2624      bValue7-=bValue0*a[7+0*BLOCKING8];
2625      bValue2-=bValue1*a[2+1*BLOCKING8];
2626      bBase[2] = bValue2;
2627      bValue3-=bValue1*a[3+1*BLOCKING8];
2628      bValue4-=bValue1*a[4+1*BLOCKING8];
2629      bValue5-=bValue1*a[5+1*BLOCKING8];
2630      bValue6-=bValue1*a[6+1*BLOCKING8];
2631      bValue7-=bValue1*a[7+1*BLOCKING8];
2632      bValue3-=bValue2*a[3+2*BLOCKING8];
2633      bBase[3] = bValue3;
2634      bValue4-=bValue2*a[4+2*BLOCKING8];
2635      bValue5-=bValue2*a[5+2*BLOCKING8];
2636      bValue6-=bValue2*a[6+2*BLOCKING8];
2637      bValue7-=bValue2*a[7+2*BLOCKING8];
2638      bValue4-=bValue3*a[4+3*BLOCKING8];
2639      bBase[4] = bValue4;
2640      bValue5-=bValue3*a[5+3*BLOCKING8];
2641      bValue6-=bValue3*a[6+3*BLOCKING8];
2642      bValue7-=bValue3*a[7+3*BLOCKING8];
2643      bValue5-=bValue4*a[5+4*BLOCKING8];
2644      bBase[5] = bValue5;
2645      bValue6-=bValue4*a[6+4*BLOCKING8];
2646      bValue7-=bValue4*a[7+4*BLOCKING8];
2647      bValue6-=bValue5*a[6+5*BLOCKING8];
2648      bBase[6] = bValue6;
2649      bValue7-=bValue5*a[7+5*BLOCKING8];
2650      bValue7-=bValue6*a[7+6*BLOCKING8];
2651      bBase[7] = bValue7;
2652#endif
2653      bBase += BLOCKING8;
2654    }
2655    bBase2 +=lda*BLOCKING8; 
2656  }
2657}
2658#define UNROLL_DTRSM 16
2659#define CILK_DTRSM 32
2660static void dtrsm0(int kkk, int first, int last,
2661                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2662{
2663  int mm=CoinMin(kkk+UNROLL_DTRSM*BLOCKING8,m);
2664  assert ((last-first)%BLOCKING8==0);
2665  if (last-first>CILK_DTRSM) {
2666    int mid=((first+last)>>4)<<3;
2667    cilk_spawn dtrsm0(kkk,first,mid,m,a,b);
2668    dtrsm0(kkk,mid,last,m,a,b);
2669    cilk_sync;
2670  } else {
2671    const long double * COIN_RESTRICT aBaseA = a+UNROLL_DTRSM*BLOCKING8X8+kkk*BLOCKING8;
2672    aBaseA += (first-mm)*BLOCKING8-BLOCKING8X8;
2673    aBaseA += m*kkk;
2674    for (int ii=first;ii<last;ii+=BLOCKING8) {
2675      aBaseA += BLOCKING8X8;
2676      const long double * COIN_RESTRICT aBaseB=aBaseA;
2677      long double * COIN_RESTRICT bBaseI = b+ii;
2678      for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
2679        long double * COIN_RESTRICT bBase = b+kk;
2680        const long double * COIN_RESTRICT aBase2 = aBaseB;//a+UNROLL_DTRSM*BLOCKING8X8+m*kk+kkk*BLOCKING8;
2681        //aBase2 += (ii-mm)*BLOCKING8;
2682        //assert (aBase2==aBaseB);
2683        aBaseB+=m*BLOCKING8;
2684#if AVX2 !=2
2685#define ALTERNATE_INNER
2686#ifndef ALTERNATE_INNER
2687        for (int k=0;k<BLOCKING8;k++) {
2688          //coin_prefetch_const(aBase2+BLOCKING8);             
2689          for (int i=0;i<BLOCKING8;i++) {
2690            bBaseI[i] -= bBase[k]*aBase2[i];
2691          }
2692          aBase2 += BLOCKING8;
2693        }
2694#else
2695        long double b0=bBaseI[0];
2696        long double b1=bBaseI[1];
2697        long double b2=bBaseI[2];
2698        long double b3=bBaseI[3];
2699        long double b4=bBaseI[4];
2700        long double b5=bBaseI[5];
2701        long double b6=bBaseI[6];
2702        long double b7=bBaseI[7];
2703        for (int k=0;k<BLOCKING8;k++) {
2704          //coin_prefetch_const(aBase2+BLOCKING8);             
2705          long double bValue=bBase[k];
2706          b0 -= bValue * aBase2[0];
2707          b1 -= bValue * aBase2[1];
2708          b2 -= bValue * aBase2[2];
2709          b3 -= bValue * aBase2[3];
2710          b4 -= bValue * aBase2[4];
2711          b5 -= bValue * aBase2[5];
2712          b6 -= bValue * aBase2[6];
2713          b7 -= bValue * aBase2[7];
2714          aBase2 += BLOCKING8;
2715        }
2716        bBaseI[0]=b0;
2717        bBaseI[1]=b1;
2718        bBaseI[2]=b2;
2719        bBaseI[3]=b3;
2720        bBaseI[4]=b4;
2721        bBaseI[5]=b5;
2722        bBaseI[6]=b6;
2723        bBaseI[7]=b7;
2724#endif
2725#else
2726        __m256d b0=_mm256_load_pd(bBaseI);
2727        __m256d b1=_mm256_load_pd(bBaseI+4);
2728        for (int k=0;k<BLOCKING8;k++) {
2729          //coin_prefetch_const(aBase2+BLOCKING8);             
2730          __m256d bb = _mm256_broadcast_sd(bBase+k);
2731          __m256d a0 = _mm256_load_pd(aBase2);
2732          b0 -= a0*bb;
2733          __m256d a1 = _mm256_load_pd(aBase2+4);
2734          b1 -= a1*bb;
2735          aBase2 += BLOCKING8;
2736        }
2737        //_mm256_store_pd ((bBaseI), (b0));
2738        *reinterpret_cast<__m256d *>(bBaseI)=b0;
2739        //_mm256_store_pd (bBaseI+4, b1);
2740        *reinterpret_cast<__m256d *>(bBaseI+4)=b1;
2741#endif
2742      }
2743    }
2744  }
2745}
2746void CoinAbcDtrsm0(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2747{
2748  assert ((m&(BLOCKING8-1))==0);
2749  // 0 Left Lower NoTranspose Unit
2750  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM*BLOCKING8) {
2751    int mm=CoinMin(m,kkk+UNROLL_DTRSM*BLOCKING8);
2752    for (int kk=kkk;kk<mm;kk+=BLOCKING8) {
2753      const long double * COIN_RESTRICT aBase2 = a+kk*(m+BLOCKING8);
2754      long double * COIN_RESTRICT bBase = b+kk;
2755      for (int k=0;k<BLOCKING8;k++) {
2756        for (int i=k+1;i<BLOCKING8;i++) {
2757          bBase[i] -= bBase[k]*aBase2[i];
2758        }
2759        aBase2 += BLOCKING8;
2760      }
2761      for (int ii=kk+BLOCKING8;ii<mm;ii+=BLOCKING8) {
2762        long double * COIN_RESTRICT bBaseI = b+ii;
2763        for (int k=0;k<BLOCKING8;k++) {
2764          //coin_prefetch_const(aBase2+BLOCKING8);             
2765          for (int i=0;i<BLOCKING8;i++) {
2766            bBaseI[i] -= bBase[k]*aBase2[i];
2767          }
2768          aBase2 += BLOCKING8;
2769        }
2770      }
2771    }
2772    dtrsm0(kkk,mm,m,m,a,b);
2773  }
2774}
2775static void dtrsm1(int kkk, int first, int last,
2776                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2777{
2778  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2779  assert ((last-first)%BLOCKING8==0);
2780  if (last-first>CILK_DTRSM) {
2781    int mid=((first+last)>>4)<<3;
2782    cilk_spawn dtrsm1(kkk,first,mid,m,a,b);
2783    dtrsm1(kkk,mid,last,m,a,b);
2784    cilk_sync;
2785  } else {
2786    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2787      long double * COIN_RESTRICT bBase2 = b+iii;
2788      const long double * COIN_RESTRICT aBaseA=
2789        a+BLOCKING8X8+BLOCKING8*iii;
2790      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2791        long double * COIN_RESTRICT bBase = b+ii;
2792        const long double * COIN_RESTRICT aBase=aBaseA+m*ii;
2793#if AVX2 !=2
2794#ifndef ALTERNATE_INNER
2795        for (int i=BLOCKING8-1;i>=0;i--) {
2796          aBase -= BLOCKING8;
2797          //coin_prefetch_const(aBase-BLOCKING8);               
2798          for (int k=BLOCKING8-1;k>=0;k--) {
2799            bBase2[k] -= bBase[i]*aBase[k];
2800          }
2801        }
2802#else
2803        long double b0=bBase2[0];
2804        long double b1=bBase2[1];
2805        long double b2=bBase2[2];
2806        long double b3=bBase2[3];
2807        long double b4=bBase2[4];
2808        long double b5=bBase2[5];
2809        long double b6=bBase2[6];
2810        long double b7=bBase2[7];
2811        for (int i=BLOCKING8-1;i>=0;i--) {
2812          aBase -= BLOCKING8;
2813          //coin_prefetch_const(aBase-BLOCKING8);               
2814          long double bValue=bBase[i];
2815          b0 -= bValue * aBase[0];
2816          b1 -= bValue * aBase[1];
2817          b2 -= bValue * aBase[2];
2818          b3 -= bValue * aBase[3];
2819          b4 -= bValue * aBase[4];
2820          b5 -= bValue * aBase[5];
2821          b6 -= bValue * aBase[6];
2822          b7 -= bValue * aBase[7];
2823        }
2824        bBase2[0]=b0;
2825        bBase2[1]=b1;
2826        bBase2[2]=b2;
2827        bBase2[3]=b3;
2828        bBase2[4]=b4;
2829        bBase2[5]=b5;
2830        bBase2[6]=b6;
2831        bBase2[7]=b7;
2832#endif
2833#else
2834        __m256d b0=_mm256_load_pd(bBase2);
2835        __m256d b1=_mm256_load_pd(bBase2+4);
2836        for (int i=BLOCKING8-1;i>=0;i--) {
2837          aBase -= BLOCKING8;
2838          //coin_prefetch_const(aBase5-BLOCKING8);             
2839          __m256d bb = _mm256_broadcast_sd(bBase+i);
2840          __m256d a0 = _mm256_load_pd(aBase);
2841          b0 -= a0*bb;
2842          __m256d a1 = _mm256_load_pd(aBase+4);
2843          b1 -= a1*bb;
2844        }
2845        //_mm256_store_pd (bBase2, b0);
2846        *reinterpret_cast<__m256d *>(bBase2)=b0;
2847        //_mm256_store_pd (bBase2+4, b1);
2848        *reinterpret_cast<__m256d *>(bBase2+4)=b1;
2849#endif
2850      }
2851    }
2852  }
2853}
2854void CoinAbcDtrsm1(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2855{
2856  assert ((m&(BLOCKING8-1))==0);
2857  // 1 Left Upper NoTranspose NonUnit
2858  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2859    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2860    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2861      const long double * COIN_RESTRICT aBase = a+m*ii+ii*BLOCKING8;
2862      long double * COIN_RESTRICT bBase = b+ii;
2863      for (int i=BLOCKING8-1;i>=0;i--) {
2864        bBase[i] *= aBase[i*(BLOCKING8+1)];
2865        for (int k=i-1;k>=0;k--) {
2866          bBase[k] -= bBase[i]*aBase[k+i*BLOCKING8];
2867        }
2868      }
2869      long double * COIN_RESTRICT bBase2 = bBase;
2870      for (int iii=ii;iii>mm;iii-=BLOCKING8) {
2871        bBase2 -= BLOCKING8;
2872        for (int i=BLOCKING8-1;i>=0;i--) {
2873          aBase -= BLOCKING8;
2874          coin_prefetch_const(aBase-BLOCKING8);         
2875          for (int k=BLOCKING8-1;k>=0;k--) {
2876            bBase2[k] -= bBase[i]*aBase[k];
2877          }
2878        }
2879      }
2880    }
2881    dtrsm1(kkk, 0, mm,  m, a, b);
2882  }
2883}
2884static void dtrsm2(int kkk, int first, int last,
2885                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2886{
2887  int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2888  assert ((last-first)%BLOCKING8==0);
2889  if (last-first>CILK_DTRSM) {
2890    int mid=((first+last)>>4)<<3;
2891    cilk_spawn dtrsm2(kkk,first,mid,m,a,b);
2892    dtrsm2(kkk,mid,last,m,a,b);
2893    cilk_sync;
2894  } else {
2895    for (int iii=last-BLOCKING8;iii>=first;iii-=BLOCKING8) {
2896      for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2897        const long double * COIN_RESTRICT bBase = b+ii;
2898        long double * COIN_RESTRICT bBase2=b+iii;
2899        const long double * COIN_RESTRICT aBase=a+ii*BLOCKING8+m*iii+BLOCKING8X8;
2900        for (int j=BLOCKING8-1;j>=0;j-=4) {
2901          long double bValue0=bBase2[j];
2902          long double bValue1=bBase2[j-1];
2903          long double bValue2=bBase2[j-2];
2904          long double bValue3=bBase2[j-3];
2905          bValue0 -= aBase[-1*BLOCKING8+7]*bBase[7];
2906          bValue1 -= aBase[-2*BLOCKING8+7]*bBase[7];
2907          bValue2 -= aBase[-3*BLOCKING8+7]*bBase[7];
2908          bValue3 -= aBase[-4*BLOCKING8+7]*bBase[7];
2909          bValue0 -= aBase[-1*BLOCKING8+6]*bBase[6];
2910          bValue1 -= aBase[-2*BLOCKING8+6]*bBase[6];
2911          bValue2 -= aBase[-3*BLOCKING8+6]*bBase[6];
2912          bValue3 -= aBase[-4*BLOCKING8+6]*bBase[6];
2913          bValue0 -= aBase[-1*BLOCKING8+5]*bBase[5];
2914          bValue1 -= aBase[-2*BLOCKING8+5]*bBase[5];
2915          bValue2 -= aBase[-3*BLOCKING8+5]*bBase[5];
2916          bValue3 -= aBase[-4*BLOCKING8+5]*bBase[5];
2917          bValue0 -= aBase[-1*BLOCKING8+4]*bBase[4];
2918          bValue1 -= aBase[-2*BLOCKING8+4]*bBase[4];
2919          bValue2 -= aBase[-3*BLOCKING8+4]*bBase[4];
2920          bValue3 -= aBase[-4*BLOCKING8+4]*bBase[4];
2921          bValue0 -= aBase[-1*BLOCKING8+3]*bBase[3];
2922          bValue1 -= aBase[-2*BLOCKING8+3]*bBase[3];
2923          bValue2 -= aBase[-3*BLOCKING8+3]*bBase[3];
2924          bValue3 -= aBase[-4*BLOCKING8+3]*bBase[3];
2925          bValue0 -= aBase[-1*BLOCKING8+2]*bBase[2];
2926          bValue1 -= aBase[-2*BLOCKING8+2]*bBase[2];
2927          bValue2 -= aBase[-3*BLOCKING8+2]*bBase[2];
2928          bValue3 -= aBase[-4*BLOCKING8+2]*bBase[2];
2929          bValue0 -= aBase[-1*BLOCKING8+1]*bBase[1];
2930          bValue1 -= aBase[-2*BLOCKING8+1]*bBase[1];
2931          bValue2 -= aBase[-3*BLOCKING8+1]*bBase[1];
2932          bValue3 -= aBase[-4*BLOCKING8+1]*bBase[1];
2933          bValue0 -= aBase[-1*BLOCKING8+0]*bBase[0];
2934          bValue1 -= aBase[-2*BLOCKING8+0]*bBase[0];
2935          bValue2 -= aBase[-3*BLOCKING8+0]*bBase[0];
2936          bValue3 -= aBase[-4*BLOCKING8+0]*bBase[0];
2937          bBase2[j]=bValue0;
2938          bBase2[j-1]=bValue1;
2939          bBase2[j-2]=bValue2;
2940          bBase2[j-3]=bValue3;
2941          aBase -= 4*BLOCKING8;
2942        }
2943      }
2944    }
2945  }
2946}
2947void CoinAbcDtrsm2(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2948{
2949  assert ((m&(BLOCKING8-1))==0);
2950  // 2 Left Lower Transpose Unit
2951  for (int kkk=m-BLOCKING8;kkk>=0;kkk-=UNROLL_DTRSM*BLOCKING8) {
2952    int mm=CoinMax(0,kkk-(UNROLL_DTRSM-1)*BLOCKING8);
2953    for (int ii=kkk;ii>=mm;ii-=BLOCKING8) {
2954      long double * COIN_RESTRICT bBase = b+ii;
2955      long double * COIN_RESTRICT bBase2 = bBase;
2956      const long double * COIN_RESTRICT aBase=a+m*ii+(ii+BLOCKING8)*BLOCKING8;
2957      for (int i=BLOCKING8-1;i>=0;i--) {
2958        aBase -= BLOCKING8;
2959        for (int k=i+1;k<BLOCKING8;k++) {
2960          bBase2[i] -= aBase[k]*bBase2[k];
2961        }
2962      }
2963      for (int iii=ii-BLOCKING8;iii>=mm;iii-=BLOCKING8) {
2964        bBase2 -= BLOCKING8;
2965        assert (bBase2==b+iii);
2966        aBase -= m*BLOCKING8;
2967        const long double * COIN_RESTRICT aBase2 = aBase+BLOCKING8X8;
2968        for (int i=BLOCKING8-1;i>=0;i--) {
2969          aBase2 -= BLOCKING8;
2970          for (int k=BLOCKING8-1;k>=0;k--) {
2971            bBase2[i] -= aBase2[k]*bBase[k];
2972          }
2973        }
2974      }
2975    }
2976    dtrsm2(kkk, 0, mm,  m, a, b);
2977  }
2978}
2979#define UNROLL_DTRSM3 16
2980#define CILK_DTRSM3 32
2981static void dtrsm3(int kkk, int first, int last,
2982                   int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
2983{
2984  //int mm=CoinMin(kkk+UNROLL_DTRSM3*BLOCKING8,m);
2985  assert ((last-first)%BLOCKING8==0);
2986  if (last-first>CILK_DTRSM3) {
2987    int mid=((first+last)>>4)<<3;
2988    cilk_spawn dtrsm3(kkk,first,mid,m,a,b);
2989    dtrsm3(kkk,mid,last,m,a,b);
2990    cilk_sync;
2991  } else {
2992    for (int kk=0;kk<kkk;kk+=BLOCKING8) {
2993      for (int ii=first;ii<last;ii+=BLOCKING8) {
2994        long double * COIN_RESTRICT bBase = b+ii;
2995        long double * COIN_RESTRICT bBase2 = b+kk;
2996        const long double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
2997        for (int j=0;j<BLOCKING8;j+=4) {
2998          long double bValue0=bBase[j];
2999          long double bValue1=bBase[j+1];
3000          long double bValue2=bBase[j+2];
3001          long double bValue3=bBase[j+3];
3002          bValue0 -= aBase[0]*bBase2[0];
3003          bValue1 -= aBase[1*BLOCKING8+0]*bBase2[0];
3004          bValue2 -= aBase[2*BLOCKING8+0]*bBase2[0];
3005          bValue3 -= aBase[3*BLOCKING8+0]*bBase2[0];
3006          bValue0 -= aBase[1]*bBase2[1];
3007          bValue1 -= aBase[1*BLOCKING8+1]*bBase2[1];
3008          bValue2 -= aBase[2*BLOCKING8+1]*bBase2[1];
3009          bValue3 -= aBase[3*BLOCKING8+1]*bBase2[1];
3010          bValue0 -= aBase[2]*bBase2[2];
3011          bValue1 -= aBase[1*BLOCKING8+2]*bBase2[2];
3012          bValue2 -= aBase[2*BLOCKING8+2]*bBase2[2];
3013          bValue3 -= aBase[3*BLOCKING8+2]*bBase2[2];
3014          bValue0 -= aBase[3]*bBase2[3];
3015          bValue1 -= aBase[1*BLOCKING8+3]*bBase2[3];
3016          bValue2 -= aBase[2*BLOCKING8+3]*bBase2[3];
3017          bValue3 -= aBase[3*BLOCKING8+3]*bBase2[3];
3018          bValue0 -= aBase[4]*bBase2[4];
3019          bValue1 -= aBase[1*BLOCKING8+4]*bBase2[4];
3020          bValue2 -= aBase[2*BLOCKING8+4]*bBase2[4];
3021          bValue3 -= aBase[3*BLOCKING8+4]*bBase2[4];
3022          bValue0 -= aBase[5]*bBase2[5];
3023          bValue1 -= aBase[1*BLOCKING8+5]*bBase2[5];
3024          bValue2 -= aBase[2*BLOCKING8+5]*bBase2[5];
3025          bValue3 -= aBase[3*BLOCKING8+5]*bBase2[5];
3026          bValue0 -= aBase[6]*bBase2[6];
3027          bValue1 -= aBase[1*BLOCKING8+6]*bBase2[6];
3028          bValue2 -= aBase[2*BLOCKING8+7]*bBase2[7];
3029          bValue3 -= aBase[3*BLOCKING8+6]*bBase2[6];
3030          bValue0 -= aBase[7]*bBase2[7];
3031          bValue1 -= aBase[1*BLOCKING8+7]*bBase2[7];
3032          bValue2 -= aBase[2*BLOCKING8+6]*bBase2[6];
3033          bValue3 -= aBase[3*BLOCKING8+7]*bBase2[7];
3034          bBase[j]=bValue0;
3035          bBase[j+1]=bValue1;
3036          bBase[j+2]=bValue2;
3037          bBase[j+3]=bValue3;
3038          aBase += 4*BLOCKING8;
3039        }
3040      }
3041    }
3042  }
3043}
3044void CoinAbcDtrsm3(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
3045{
3046  assert ((m&(BLOCKING8-1))==0);
3047  // 3 Left Upper Transpose NonUnit
3048  for (int kkk=0;kkk<m;kkk+=UNROLL_DTRSM3*BLOCKING8) {
3049    int mm=CoinMin(m,kkk+UNROLL_DTRSM3*BLOCKING8);
3050    dtrsm3(kkk, kkk, mm,  m, a, b);
3051    for (int ii=kkk;ii<mm;ii+=BLOCKING8) {
3052      long double * COIN_RESTRICT bBase = b+ii;
3053      for (int kk=kkk;kk<ii;kk+=BLOCKING8) {
3054        long double * COIN_RESTRICT bBase2 = b+kk;
3055        const long double * COIN_RESTRICT aBase=a+ii*m+kk*BLOCKING8;
3056        for (int i=0;i<BLOCKING8;i++) {
3057          for (int k=0;k<BLOCKING8;k++) {
3058            bBase[i] -= aBase[k]*bBase2[k];
3059          }
3060          aBase += BLOCKING8;
3061        }
3062      }
3063      const long double * COIN_RESTRICT aBase=a+ii*m+ii*BLOCKING8;
3064      for (int i=0;i<BLOCKING8;i++) {
3065        for (int k=0;k<i;k++) {
3066          bBase[i] -= aBase[k]*bBase[k];
3067        }
3068        bBase[i] = bBase[i]*aBase[i];
3069        aBase += BLOCKING8;
3070      }
3071    }
3072  }
3073}
3074static void 
3075CoinAbcDlaswp(int n, long double * COIN_RESTRICT a,int lda,int start,int end, int * ipiv)
3076{
3077  assert ((n&(BLOCKING8-1))==0);
3078  long double * COIN_RESTRICT aThis3 = a;
3079  for (int j=0;j<n;j+=BLOCKING8) {
3080    for (int i=start;i<end;i++) {
3081      int ip=ipiv[i];
3082      if (ip!=i) {
3083        long double * COIN_RESTRICT aTop=aThis3+j*lda;
3084        int iBias = ip/BLOCKING8;
3085        ip-=iBias*BLOCKING8;
3086        long double * COIN_RESTRICT aNotTop=aTop+iBias*BLOCKING8X8;
3087        iBias = i/BLOCKING8;
3088        int i2=i-iBias*BLOCKING8;
3089        aTop += iBias*BLOCKING8X8;
3090        for (int k=0;k<BLOCKING8;k++) {
3091          long double temp=aTop[i2];
3092          aTop[i2]=aNotTop[ip];
3093          aNotTop[ip]=temp;
3094          aTop += BLOCKING8;
3095          aNotTop += BLOCKING8;
3096        }
3097      }
3098    }
3099  }
3100}
3101extern void CoinAbcDgemm(int m, int n, int k, long double * COIN_RESTRICT a,int lda,
3102                          long double * COIN_RESTRICT b,long double * COIN_RESTRICT c
3103#if ABC_PARALLEL==2
3104                          ,int parallelMode
3105#endif
3106                         );
3107static int CoinAbcDgetrf2(int m, int n, long double * COIN_RESTRICT a, int * ipiv)
3108{
3109  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
3110  assert (n==BLOCKING8);
3111  int dimension=CoinMin(m,n);
3112  /* entry for column j and row i (when multiple of BLOCKING8)
3113     is at aBlocked+j*m+i*BLOCKING8
3114  */
3115  assert (dimension==BLOCKING8);
3116  //printf("Dgetrf2 m=%d n=%d\n",m,n);
3117  long double * COIN_RESTRICT aThis3 = a;
3118  for (int j=0;j<dimension;j++) {
3119    int pivotRow=-1;
3120    long double largest=0.0;
3121    long double * COIN_RESTRICT aThis2 = aThis3+j*BLOCKING8;
3122    // this block
3123    int pivotRow2=-1;
3124    for (int i=j;i<BLOCKING8;i++) {
3125      long double value=fabsl(aThis2[i]);
3126      if (value>largest) {
3127        largest=value;
3128        pivotRow2=i;
3129      }
3130    }
3131    // other blocks
3132    int iBias=0;
3133    for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
3134      aThis2+=BLOCKING8*BLOCKING8;
3135      for (int i=0;i<BLOCKING8;i++) {
3136        long double value=fabsl(aThis2[i]);
3137        if (value>largest) {
3138          largest=value;
3139          pivotRow2=i;
3140          iBias=ii;
3141        }
3142      }
3143    }
3144    pivotRow=pivotRow2+iBias;
3145    ipiv[j]=pivotRow;
3146    if (largest) {
3147      if (j!=pivotRow) {
3148        long double * COIN_RESTRICT aTop=aThis3;
3149        long double * COIN_RESTRICT aNotTop=aThis3+iBias*BLOCKING8;
3150        for (int i=0;i<n;i++) {
3151          long double value = aTop[j];
3152          aTop[j]=aNotTop[pivotRow2];
3153          aNotTop[pivotRow2]=value;
3154          aTop += BLOCKING8;
3155          aNotTop += BLOCKING8;
3156        }
3157      }
3158      aThis2 = aThis3+j*BLOCKING8;
3159      long double pivotMultiplier = 1.0 / aThis2[j];
3160      aThis2[j] = pivotMultiplier;
3161      // Do L
3162      // this block
3163      for (int i=j+1;i<BLOCKING8;i++) 
3164        aThis2[i] *= pivotMultiplier;
3165      // other blocks
3166      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
3167        aThis2+=BLOCKING8*BLOCKING8;
3168        for (int i=0;i<BLOCKING8;i++) 
3169          aThis2[i] *= pivotMultiplier;
3170      }
3171      // update rest
3172      long double * COIN_RESTRICT aPut;
3173      aThis2 = aThis3+j*BLOCKING8;
3174      aPut = aThis2+BLOCKING8;
3175      long double * COIN_RESTRICT aPut2 = aPut;
3176      // this block
3177      for (int i=j+1;i<BLOCKING8;i++) {
3178        long double value = aPut[j];
3179        for (int k=j+1;k<BLOCKING8;k++) 
3180          aPut[k] -= value*aThis2[k];
3181        aPut += BLOCKING8;
3182      }
3183      // other blocks
3184      for (int ii=BLOCKING8;ii<m;ii+=BLOCKING8) {
3185        aThis2 += BLOCKING8*BLOCKING8;
3186        aPut = aThis2+BLOCKING8;
3187        long double * COIN_RESTRICT aPut2a = aPut2;
3188        for (int i=j+1;i<BLOCKING8;i++) {
3189          long double value = aPut2a[j];
3190          for (int k=0;k<BLOCKING8;k++) 
3191            aPut[k] -= value*aThis2[k];
3192          aPut += BLOCKING8;
3193          aPut2a += BLOCKING8;
3194        }
3195      }
3196    } else {
3197      return 1;
3198    }
3199  }
3200  return 0;
3201}
3202
3203int 
3204CoinAbcDgetrf(int m, int n, long double * COIN_RESTRICT a, int lda, int * ipiv
3205#if ABC_PARALLEL==2
3206                          ,int parallelMode
3207#endif
3208)
3209{
3210  assert (m==n);
3211  assert ((m&(BLOCKING8-1))==0&&(n&(BLOCKING8-1))==0);
3212  if (m<BLOCKING8) {
3213    return CoinAbcDgetrf2(m,n,a,ipiv);
3214  } else {
3215    for (int j=0;j<n;j+=BLOCKING8) {
3216      int start=j;
3217      int newSize=CoinMin(BLOCKING8,n-j);
3218      int end=j+newSize;
3219      int returnCode=CoinAbcDgetrf2(m-start,newSize,a+(start*lda+start*BLOCKING8),
3220                                     ipiv+start);
3221      if (!returnCode) {
3222        // adjust
3223        for (int k=start;k<end;k++)
3224          ipiv[k]+=start;
3225        // swap 0<start
3226        CoinAbcDlaswp(start,a,lda,start,end,ipiv);
3227        if (end<n) {
3228          // swap >=end
3229          CoinAbcDlaswp(n-end,a+end*lda,lda,start,end,ipiv);
3230          CoinAbcDtrsmFactor(newSize,n-end,a+(start*lda+start*BLOCKING8),lda);
3231          CoinAbcDgemm(n-end,n-end,newSize,
3232                        a+start*lda+end*BLOCKING8,lda,
3233                        a+end*lda+start*BLOCKING8,a+end*lda+end*BLOCKING8
3234#if ABC_PARALLEL==2
3235                          ,parallelMode
3236#endif
3237                        );
3238        }
3239      } else {
3240        return returnCode;
3241      }
3242    }
3243  }
3244  return 0;
3245}
3246void
3247CoinAbcDgetrs(char trans,int m, long double * COIN_RESTRICT a, long double * COIN_RESTRICT work)
3248{
3249  assert ((m&(BLOCKING8-1))==0);
3250  if (trans=='N') {
3251    //CoinAbcDlaswp1(work,m,ipiv);
3252    CoinAbcDtrsm0(m,a,work);
3253    CoinAbcDtrsm1(m,a,work);
3254  } else {
3255    assert (trans=='T');
3256    CoinAbcDtrsm3(m,a,work);
3257    CoinAbcDtrsm2(m,a,work);
3258    //CoinAbcDlaswp1Backwards(work,m,ipiv);
3259  }
3260}
3261#endif
3262#endif
Note: See TracBrowser for help on using the repository browser.