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

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