Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/AbcSimplexParallel.cpp

    r2070 r2385  
    2323#include <iostream>
    2424//#undef AVX2
    25 #define _mm256_broadcast_sd(x) static_cast<__m256d> (__builtin_ia32_vbroadcastsd256 (x))
     25#define _mm256_broadcast_sd(x) static_cast< __m256d >(__builtin_ia32_vbroadcastsd256(x))
    2626#define _mm256_load_pd(x) *(__m256d *)(x)
    27 #define _mm256_store_pd (s, x)  *((__m256d *)s)=x
     27#define _mm256_store_pd (s, x) * ((__m256d *)s) = x
    2828//#define ABC_DEBUG 2
    2929/* Reasons to come out:
     
    3535   +3 max iterations
    3636*/
    37 int
    38 AbcSimplexDual::whileIteratingSerial()
     37int AbcSimplexDual::whileIteratingSerial()
    3938{
    4039  /* arrays
     
    5655  int returnCode = -1;
    5756#define DELAYED_UPDATE
    58   arrayForBtran_=0;
     57  arrayForBtran_ = 0;
    5958  setUsedArray(arrayForBtran_);
    60   arrayForFtran_=1;
     59  arrayForFtran_ = 1;
    6160  setUsedArray(arrayForFtran_);
    62   arrayForFlipBounds_=2;
     61  arrayForFlipBounds_ = 2;
    6362  setUsedArray(arrayForFlipBounds_);
    64   arrayForTableauRow_=3;
     63  arrayForTableauRow_ = 3;
    6564  setUsedArray(arrayForTableauRow_);
    66   arrayForDualColumn_=4;
     65  arrayForDualColumn_ = 4;
    6766  setUsedArray(arrayForDualColumn_);
    68   arrayForReplaceColumn_=4; //4;
    69   arrayForFlipRhs_=5;
     67  arrayForReplaceColumn_ = 4; //4;
     68  arrayForFlipRhs_ = 5;
    7069  setUsedArray(arrayForFlipRhs_);
    7170  dualPivotRow();
    72   lastPivotRow_=pivotRow_;
     71  lastPivotRow_ = pivotRow_;
    7372  if (pivotRow_ >= 0) {
    7473    // we found a pivot row
     
    8483        2 - don't take and break
    8584      */
    86       stateOfIteration_=0;
    87       returnCode=-1;
     85      stateOfIteration_ = 0;
     86      returnCode = -1;
    8887      // put row of tableau in usefulArray[arrayForTableauRow_]
    8988      /*
     
    9493      */
    9594      //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;
     95      double saveAcceptable = acceptablePivot_;
     96      if (largestPrimalError_ > 1.0e-5 || largestDualError_ > 1.0e-5) {
     97        if (!abcFactorization_->pivots())
     98          acceptablePivot_ *= 1.0e2;
     99        else if (abcFactorization_->pivots() < 5)
     100          acceptablePivot_ *= 1.0e1;
    102101      }
    103102      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         }
     103      acceptablePivot_ = saveAcceptable;
     104      if ((stateOfProblem_ & VALUES_PASS) != 0) {
     105        // see if can just move dual
     106        if (fabs(upperTheta_ - fabs(abcDj_[sequenceOut_])) < dualTolerance_) {
     107          stateOfIteration_ = -1;
     108        }
    110109      }
    111110      //usefulArray_[arrayForTableauRow_].sortPacked();
     
    115114      //usefulArray_[arrayForDualColumn].print();
    116115      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         }
     116        // get sequenceIn_
     117        dualPivotColumn();
     118        if (sequenceIn_ >= 0) {
     119          // normal iteration
     120          // update the incoming column (and do weights)
     121          getTableauColumnFlipAndStartReplaceSerial();
     122        } else if (stateOfIteration_ != -1) {
     123          stateOfIteration_ = 2;
     124        }
    126125      }
    127126      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_);
     127        //      assert (stateOfIteration_!=1);
     128        int whatNext = 0;
     129        whatNext = housekeeping();
     130        if (whatNext == 1) {
     131          problemStatus_ = -2; // refactorize
     132        } else if (whatNext == 2) {
     133          // maximum iterations or equivalent
     134          problemStatus_ = 3;
     135          returnCode = 3;
     136          stateOfIteration_ = 2;
     137        }
     138        if (problemStatus_ == -1) {
     139          replaceColumnPart3();
     140          //clearArrays(arrayForReplaceColumn_);
    142141#if ABC_DEBUG
    143           checkArrays();
    144 #endif
    145             updateDualsInDual();
    146             abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
    147                                                              movement_);
     142          checkArrays();
     143#endif
     144          updateDualsInDual();
     145          abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_],
     146            movement_);
    148147#if ABC_DEBUG
    149           checkArrays();
    150 #endif
    151         } else {
    152           abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
    153           //clearArrays(arrayForBtran_);
    154           //clearArrays(arrayForFtran_);
    155         }
     148          checkArrays();
     149#endif
     150        } else {
     151          abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_]);
     152          //clearArrays(arrayForBtran_);
     153          //clearArrays(arrayForFtran_);
     154        }
    156155      } 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;
     156        if (stateOfIteration_ < 0) {
     157          // move dual in dual values pass
     158          theta_ = abcDj_[sequenceOut_];
     159          updateDualsInDual();
     160          abcDj_[sequenceOut_] = 0.0;
     161          sequenceOut_ = -1;
     162        }
     163        // clear all arrays
     164        clearArrays(-1);
     165        //sequenceIn_=-1;
     166        //sequenceOut_=-1;
    168167      }
    169168      // Check event
    170169      {
    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         }
     170        int status = eventHandler_->event(ClpEventHandler::endOfIteration);
     171        if (status >= 0) {
     172          problemStatus_ = 5;
     173          secondaryStatus_ = ClpEventHandler::endOfIteration;
     174          returnCode = 4;
     175          stateOfIteration_ = 2;
     176        }
    178177      }
    179178      // 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;
     179      if (sequenceIn_ < 0 && sequenceOut_ >= 0) {
     180        // no incoming column is valid
     181        returnCode = noPivotColumn();
     182      }
     183      if (stateOfIteration_ == 2) {
     184        sequenceIn_ = -1;
     185        break;
    187186      }
    188187      swapPrimalStuff();
    189       if (problemStatus_!=-1) {
    190         break;
     188      if (problemStatus_ != -1) {
     189        break;
    191190      }
    192191      // dualRow will go to virtual row pivot choice algorithm
    193192      // make sure values pass off if it should be
    194193      // use Btran array and clear inside dualPivotRow (if used)
    195       int lastSequenceOut=sequenceOut_;
    196       int lastDirectionOut=directionOut_;
     194      int lastSequenceOut = sequenceOut_;
     195      int lastDirectionOut = directionOut_;
    197196      dualPivotRow();
    198       lastPivotRow_=pivotRow_;
     197      lastPivotRow_ = pivotRow_;
    199198      if (pivotRow_ >= 0) {
    200         // we found a pivot row
    201         // create as packed
    202         createDualPricingVectorSerial();
    203         swapDualStuff(lastSequenceOut,lastDirectionOut);
    204         // next pivot
     199        // we found a pivot row
     200        // create as packed
     201        createDualPricingVectorSerial();
     202        swapDualStuff(lastSequenceOut, lastDirectionOut);
     203        // next pivot
    205204      } else {
    206         // no pivot row
    207         clearArrays(-1);
    208         returnCode=noPivotRow();
    209         break;
     205        // no pivot row
     206        clearArrays(-1);
     207        returnCode = noPivotRow();
     208        break;
    210209      }
    211210    } while (problemStatus_ == -1);
     
    217216    // no pivot row on first try
    218217    clearArrays(-1);
    219     returnCode=noPivotRow();
     218    returnCode = noPivotRow();
    220219  }
    221220  //printStuff();
     
    226225}
    227226// Create dual pricing vector
    228 void
    229 AbcSimplexDual::createDualPricingVectorSerial()
     227void AbcSimplexDual::createDualPricingVectorSerial()
    230228{
    231229#ifndef NDEBUG
    232 #if ABC_NORMAL_DEBUG>3
    233     checkArrays();
     230#if ABC_NORMAL_DEBUG > 3
     231  checkArrays();
    234232#endif
    235233#endif
     
    248246  sequenceIn_ = -1;
    249247}
    250 void
    251 AbcSimplexDual::getTableauColumnPart1Serial()
     248void AbcSimplexDual::getTableauColumnPart1Serial()
    252249{
    253250#if ABC_DEBUG
    254251  {
    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;
     252    const double *work = usefulArray_[arrayForTableauRow_].denseVector();
     253    int number = usefulArray_[arrayForTableauRow_].getNumElements();
     254    const int *which = usefulArray_[arrayForTableauRow_].getIndices();
     255    for (int i = 0; i < number; i++) {
     256      if (which[i] == sequenceIn_) {
     257        assert(alpha_ == work[i]);
     258        break;
    262259      }
    263260    }
     
    269266  // Array[0] may be needed until updateWeights2
    270267  // and update dual weights (can do in parallel - with extra array)
    271   alpha_ = abcDualRowPivot_->updateWeights1(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
    272 }
    273 int
    274 AbcSimplexDual::getTableauColumnFlipAndStartReplaceSerial()
     268  alpha_ = abcDualRowPivot_->updateWeights1(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_]);
     269}
     270int AbcSimplexDual::getTableauColumnFlipAndStartReplaceSerial()
    275271{
    276272  // move checking stuff down into called functions
    277273  int numberFlipped;
    278   numberFlipped=flipBounds();
     274  numberFlipped = flipBounds();
    279275  // update the incoming column
    280276  getTableauColumnPart1Serial();
     
    285281  //usefulArray_[arrayForTableauRow_].compact();
    286282  // returns 3 if skip this iteration and re-factorize
    287 /*
     283  /*
    288284  return code
    289285  0 - OK
     
    292288  4 - break and say infeasible
    293289 */
    294   int returnCode=0;
     290  int returnCode = 0;
    295291  // amount primal will move
    296292  movement_ = -dualOut_ * directionOut_ / alpha_;
    297293  // see if update stable
    298 #if ABC_NORMAL_DEBUG>3
     294#if ABC_NORMAL_DEBUG > 3
    299295  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_)) {
     296    double ft = ftAlpha_ * abcFactorization_->pivotRegion()[pivotRow_];
     297    double diff1 = fabs(alpha_ - btranAlpha_);
     298    double diff2 = fabs(fabs(alpha_) - fabs(ft));
     299    double diff3 = fabs(fabs(ft) - fabs(btranAlpha_));
     300    double largest = CoinMax(CoinMax(diff1, diff2), diff3);
     301    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n",
     302      btranAlpha_, alpha_, ft, largest);
     303    if (largest > 0.001 * fabs(alpha_)) {
    308304      printf("bad\n");
    309305    }
    310306  }
    311307#endif
    312   int numberPivots=abcFactorization_->pivots();
     308  int numberPivots = abcFactorization_->pivots();
    313309  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
    314310  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
     
    316312  if (largestPrimalError_ > 10.0)
    317313    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_))) {
     314  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > checkValue * (1.0 + fabs(alpha_))) {
    320315    handler_->message(CLP_DUAL_CHECK, messages_)
    321316      << btranAlpha_
     
    327322      test = 1.0e-1 * fabs(alpha_);
    328323    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) {
     324      test = 10.0 * checkValue; //1.0e-4 * (1.0 + fabs(alpha_));
     325    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > test);
     326    double derror = CoinMin(fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)), 1.0) * 0.9999e7;
     327    int error = 0;
     328    while (derror > 1.0) {
    335329      error++;
    336330      derror *= 0.1;
    337331    }
    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_)
     332    int frequency[8] = { 99999999, 100, 10, 2, 1, 1, 1, 1 };
     333    int newFactorFrequency = CoinMin(forceFactorization_, frequency[error]);
     334#if ABC_NORMAL_DEBUG > 0
     335    if (newFactorFrequency < forceFactorization_)
    342336      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;
     337        fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)), numberPivots,
     338        forceFactorization_, newFactorFrequency);
     339#endif
     340    if (!numberPivots && fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)) > 0.5e-6
     341      && abcFactorization_->pivotTolerance() < 0.5)
     342      abcFactorization_->saferTolerances(1.0, -1.03);
     343    forceFactorization_ = newFactorFrequency;
    350344
    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_)));
     345#if ABC_NORMAL_DEBUG > 0
     346    if (fabs(btranAlpha_ + alpha_) < 1.0e-5 * (1.0 + fabs(alpha_))) {
     347      printf("diff (%g,%g) %g check %g\n", btranAlpha_, alpha_, fabs(btranAlpha_ - alpha_), checkValue * (1.0 + fabs(alpha_)));
    355348    }
    356349#endif
    357350    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
     351      if (needFlag && numberPivots < 10) {
     352        // need to reject something
     353        assert(sequenceOut_ >= 0);
     354        char x = isColumn(sequenceOut_) ? 'C' : 'R';
     355        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     356          << x << sequenceWithin(sequenceOut_)
     357          << CoinMessageEol;
     358#if ABC_NORMAL_DEBUG > 0
     359        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
     360          btranAlpha_, alpha_, numberPivots);
     361#endif
     362        // Make safer?
     363        abcFactorization_->saferTolerances(1.0, -1.03);
     364        setFlagged(sequenceOut_);
     365        // so can't happen immediately again
     366        sequenceOut_ = -1;
     367        lastBadIteration_ = numberIterations_; // say be more cautious
    375368      }
    376369      // Make safer?
     
    379372      problemStatus_ = -2; // factorize now
    380373      returnCode = -2;
    381       stateOfIteration_=2;
     374      stateOfIteration_ = 2;
    382375    } else {
    383376      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         }
     377        assert(sequenceOut_ >= 0);
     378        // need to reject something
     379        char x = isColumn(sequenceOut_) ? 'C' : 'R';
     380        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     381          << x << sequenceWithin(sequenceOut_)
     382          << CoinMessageEol;
     383#if ABC_NORMAL_DEBUG > 3
     384        printf("flag a %g %g\n", btranAlpha_, alpha_);
     385#endif
     386        setFlagged(sequenceOut_);
     387        // so can't happen immediately again
     388        sequenceOut_ = -1;
     389        //abcProgress_.clearBadTimes();
     390        lastBadIteration_ = numberIterations_; // say be more cautious
     391        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
     392          //printf("I think should declare infeasible\n");
     393          problemStatus_ = 1;
     394          returnCode = 1;
     395          stateOfIteration_ = 2;
     396        } else {
     397          stateOfIteration_ = 1;
     398        }
     399        // Make safer?
     400        if (abcFactorization_->pivotTolerance() < 0.999 && stateOfIteration_ == 1) {
     401          // change tolerance and re-invert
     402          abcFactorization_->saferTolerances(1.0, -1.03);
     403          problemStatus_ = -6; // factorize now
     404          returnCode = -2;
     405          stateOfIteration_ = 2;
     406        }
    414407      }
    415408    }
     
    417410  if (!stateOfIteration_) {
    418411    // check update
    419     //int updateStatus = 
    420 /*
     412    //int updateStatus =
     413    /*
    421414  returns
    422415  0 - OK
     
    426419  5 - take iteration then re-factorize because of memory
    427420 */
    428     int status=checkReplace();
    429     if (status&&!returnCode)
    430       returnCode=status;
     421    int status = checkReplace();
     422    if (status && !returnCode)
     423      returnCode = status;
    431424  }
    432425  //clearArrays(arrayForFlipRhs_);
    433   if (stateOfIteration_&&numberFlipped) {
     426  if (stateOfIteration_ && numberFlipped) {
    434427    //usefulArray_[arrayForTableauRow_].compact();
    435428    // move variables back
     
    439432  return returnCode;
    440433}
    441 #if ABC_PARALLEL==1
     434#if ABC_PARALLEL == 1
    442435/* Reasons to come out:
    443436   -1 iterations etc
     
    448441   +3 max iterations
    449442*/
    450 int
    451 AbcSimplexDual::whileIteratingThread()
     443int AbcSimplexDual::whileIteratingThread()
    452444{
    453445  /* arrays
     
    469461  int returnCode = -1;
    470462#define DELAYED_UPDATE
    471   arrayForBtran_=0;
     463  arrayForBtran_ = 0;
    472464  setUsedArray(arrayForBtran_);
    473   arrayForFtran_=1;
     465  arrayForFtran_ = 1;
    474466  setUsedArray(arrayForFtran_);
    475   arrayForFlipBounds_=2;
     467  arrayForFlipBounds_ = 2;
    476468  setUsedArray(arrayForFlipBounds_);
    477   arrayForTableauRow_=3;
     469  arrayForTableauRow_ = 3;
    478470  setUsedArray(arrayForTableauRow_);
    479   arrayForDualColumn_=4;
     471  arrayForDualColumn_ = 4;
    480472  setUsedArray(arrayForDualColumn_);
    481   arrayForReplaceColumn_=4; //4;
    482   arrayForFlipRhs_=5;
     473  arrayForReplaceColumn_ = 4; //4;
     474  arrayForFlipRhs_ = 5;
    483475  setUsedArray(arrayForFlipRhs_);
    484476  dualPivotRow();
    485   lastPivotRow_=pivotRow_;
     477  lastPivotRow_ = pivotRow_;
    486478  if (pivotRow_ >= 0) {
    487479    // we found a pivot row
     
    497489        2 - don't take and break
    498490      */
    499       stateOfIteration_=0;
    500       returnCode=-1;
     491      stateOfIteration_ = 0;
     492      returnCode = -1;
    501493      // put row of tableau in usefulArray[arrayForTableauRow_]
    502494      /*
     
    507499      */
    508500      //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;
     501      double saveAcceptable = acceptablePivot_;
     502      if (largestPrimalError_ > 1.0e-5 || largestDualError_ > 1.0e-5) {
     503        if (!abcFactorization_->pivots())
     504          acceptablePivot_ *= 1.0e2;
     505        else if (abcFactorization_->pivots() < 5)
     506          acceptablePivot_ *= 1.0e1;
    515507      }
    516508      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         }
     509      acceptablePivot_ = saveAcceptable;
     510      if ((stateOfProblem_ & VALUES_PASS) != 0) {
     511        // see if can just move dual
     512        if (fabs(upperTheta_ - fabs(abcDj_[sequenceOut_])) < dualTolerance_) {
     513          stateOfIteration_ = -1;
     514        }
    523515      }
    524516      //usefulArray_[arrayForTableauRow_].sortPacked();
     
    527519      //usefulArray_[arrayForDualColumn].sortPacked();
    528520      //usefulArray_[arrayForDualColumn].print();
    529       if (parallelMode_!=0) {
    530         stopStart_=1+32*1; // Just first thread for updateweights
    531         startParallelStuff(2);
     521      if (parallelMode_ != 0) {
     522        stopStart_ = 1 + 32 * 1; // Just first thread for updateweights
     523        startParallelStuff(2);
    532524      }
    533525      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);
     526        // get sequenceIn_
     527        dualPivotColumn();
     528        if (sequenceIn_ >= 0) {
     529          // normal iteration
     530          // update the incoming column (and do weights if serial)
     531          getTableauColumnFlipAndStartReplaceThread();
     532          //usleep(1000);
     533        } else if (stateOfIteration_ != -1) {
     534          stateOfIteration_ = 2;
     535        }
     536      }
     537      if (parallelMode_ != 0) {
     538        // do sync here
     539        stopParallelStuff(2);
    548540      }
    549541      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         }
     542        //      assert (stateOfIteration_!=1);
     543        int whatNext = 0;
     544        whatNext = housekeeping();
     545        if (whatNext == 1) {
     546          problemStatus_ = -2; // refactorize
     547        } else if (whatNext == 2) {
     548          // maximum iterations or equivalent
     549          problemStatus_ = 3;
     550          returnCode = 3;
     551          stateOfIteration_ = 2;
     552        }
    561553      } 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);
     554        if (stateOfIteration_ < 0) {
     555          // move dual in dual values pass
     556          theta_ = abcDj_[sequenceOut_];
     557          updateDualsInDual();
     558          abcDj_[sequenceOut_] = 0.0;
     559          sequenceOut_ = -1;
     560        }
     561        // clear all arrays
     562        clearArrays(-1);
    571563      }
    572564      // 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;
     565      if (sequenceIn_ < 0 && sequenceOut_ >= 0) {
     566        // no incoming column is valid
     567        returnCode = noPivotColumn();
     568      }
     569      if (stateOfIteration_ == 2) {
     570        sequenceIn_ = -1;
     571        break;
     572      }
     573      if (problemStatus_ != -1) {
     574        abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_]);
     575        swapPrimalStuff();
     576        break;
    585577      }
    586578#if ABC_DEBUG
    587579      checkArrays();
    588580#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         }
     581      if (stateOfIteration_ != -1) {
     582        assert(stateOfIteration_ == 0); // if 1 why are we here
     583        // can do these in parallel
     584        if (parallelMode_ == 0) {
     585          replaceColumnPart3();
     586          updateDualsInDual();
     587          abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_],
     588            movement_);
     589        } else {
     590          stopStart_ = 1 + 32 * 1;
     591          startParallelStuff(3);
     592          if (parallelMode_ > 1) {
     593            stopStart_ = 2 + 32 * 2;
     594            startParallelStuff(9);
     595          } else {
     596            replaceColumnPart3();
     597          }
     598          abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_],
     599            movement_);
     600        }
    609601      }
    610602#if ABC_DEBUG
     
    615607      // make sure values pass off if it should be
    616608      // use Btran array and clear inside dualPivotRow (if used)
    617       int lastSequenceOut=sequenceOut_;
    618       int lastDirectionOut=directionOut_;
     609      int lastSequenceOut = sequenceOut_;
     610      int lastDirectionOut = directionOut_;
    619611      // 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);
     612      if (parallelMode_ > 1) {
     613        stopStart_ = 2 + 32 * 2; // stop factorization update
     614        //stopStart_=3+32*3; // stop all
     615        stopParallelStuff(9);
    624616      }
    625617      dualPivotRow();
    626       lastPivotRow_=pivotRow_;
     618      lastPivotRow_ = pivotRow_;
    627619      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         }
     620        // we found a pivot row
     621        // create as packed
     622        createDualPricingVectorThread();
     623        swapDualStuff(lastSequenceOut, lastDirectionOut);
     624        // next pivot
     625        // redo to test on parallelMode_
     626        if (parallelMode_ != 0) {
     627          stopStart_ = 1 + 32 * 1; // stop dual update
     628          stopParallelStuff(3);
     629        }
    638630      } 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;
     631        // no pivot row
     632        // redo to test on parallelMode_
     633        if (parallelMode_ != 0) {
     634          stopStart_ = 1 + 32 * 1; // stop dual update
     635          stopParallelStuff(3);
     636        }
     637        clearArrays(-1);
     638        returnCode = noPivotRow();
     639        break;
    648640      }
    649641    } while (problemStatus_ == -1);
     
    655647    // no pivot row on first try
    656648    clearArrays(-1);
    657     returnCode=noPivotRow();
     649    returnCode = noPivotRow();
    658650  }
    659651  //printStuff();
     
    664656}
    665657// Create dual pricing vector
    666 void
    667 AbcSimplexDual::createDualPricingVectorThread()
     658void AbcSimplexDual::createDualPricingVectorThread()
    668659{
    669660#ifndef NDEBUG
    670 #if ABC_NORMAL_DEBUG>3
    671     checkArrays();
     661#if ABC_NORMAL_DEBUG > 3
     662  checkArrays();
    672663#endif
    673664#endif
     
    686677  sequenceIn_ = -1;
    687678}
    688 void
    689 AbcSimplexDual::getTableauColumnPart1Thread()
     679void AbcSimplexDual::getTableauColumnPart1Thread()
    690680{
    691681#if ABC_DEBUG
    692682  {
    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;
     683    const double *work = usefulArray_[arrayForTableauRow_].denseVector();
     684    int number = usefulArray_[arrayForTableauRow_].getNumElements();
     685    const int *which = usefulArray_[arrayForTableauRow_].getIndices();
     686    for (int i = 0; i < number; i++) {
     687      if (which[i] == sequenceIn_) {
     688        assert(alpha_ == work[i]);
     689        break;
    700690      }
    701691    }
     
    707697  // Array[0] may be needed until updateWeights2
    708698  // and update dual weights (can do in parallel - with extra array)
    709   if (parallelMode_!=0) {
     699  if (parallelMode_ != 0) {
    710700    abcFactorization_->updateColumnFTPart1(usefulArray_[arrayForFtran_]);
    711701    // pivot element
    712702    //alpha_ = usefulArray_[arrayForFtran_].denseVector()[pivotRow_];
    713     } else {
    714     alpha_ = abcDualRowPivot_->updateWeights1(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
    715   }
    716 }
    717 int
    718 AbcSimplexDual::getTableauColumnFlipAndStartReplaceThread()
     703  } else {
     704    alpha_ = abcDualRowPivot_->updateWeights1(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_]);
     705  }
     706}
     707int AbcSimplexDual::getTableauColumnFlipAndStartReplaceThread()
    719708{
    720709  // move checking stuff down into called functions
     
    722711  int numberFlipped;
    723712
    724   if (parallelMode_==0) {
    725     numberFlipped=flipBounds();
     713  if (parallelMode_ == 0) {
     714    numberFlipped = flipBounds();
    726715    // update the incoming column
    727716    getTableauColumnPart1Thread();
     
    732721  } else {
    733722    // save stopStart
    734     int saveStopStart=stopStart_;
    735     if (parallelMode_>1) {
     723    int saveStopStart = stopStart_;
     724    if (parallelMode_ > 1) {
    736725      // we can flip immediately
    737       stopStart_=2+32*2; // just thread 1
     726      stopStart_ = 2 + 32 * 2; // just thread 1
    738727      startParallelStuff(8);
    739728      // update the incoming column
    740729      getTableauColumnPart1Thread();
    741       stopStart_=4+32*4; // just thread 2
     730      stopStart_ = 4 + 32 * 4; // just thread 2
    742731      startParallelStuff(7);
    743732      getTableauColumnPart2();
    744       stopStart_=6+32*6; // just threads 1 and 2
     733      stopStart_ = 6 + 32 * 6; // just threads 1 and 2
    745734      stopParallelStuff(8);
    746735      //ftAlpha_=threadInfo_[2].result;
    747736    } 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;
     737      // just one extra thread - do replace
     738      // update the incoming column
     739      getTableauColumnPart1Thread();
     740      stopStart_ = 2 + 32 * 2; // just thread 1
     741      startParallelStuff(7);
     742      getTableauColumnPart2();
     743      numberFlipped = flipBounds();
     744      stopParallelStuff(8);
     745      //ftAlpha_=threadInfo_[1].result;
     746    }
     747    stopStart_ = saveStopStart;
    759748  }
    760749  //usefulArray_[arrayForTableauRow_].compact();
    761750  // returns 3 if skip this iteration and re-factorize
    762 /*
     751  /*
    763752  return code
    764753  0 - OK
     
    767756  4 - break and say infeasible
    768757 */
    769   int returnCode=0;
     758  int returnCode = 0;
    770759  // amount primal will move
    771760  movement_ = -dualOut_ * directionOut_ / alpha_;
    772761  // see if update stable
    773 #if ABC_NORMAL_DEBUG>3
     762#if ABC_NORMAL_DEBUG > 3
    774763  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_)) {
     764    double ft = ftAlpha_ * abcFactorization_->pivotRegion()[pivotRow_];
     765    double diff1 = fabs(alpha_ - btranAlpha_);
     766    double diff2 = fabs(fabs(alpha_) - fabs(ft));
     767    double diff3 = fabs(fabs(ft) - fabs(btranAlpha_));
     768    double largest = CoinMax(CoinMax(diff1, diff2), diff3);
     769    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n",
     770      btranAlpha_, alpha_, ft, largest);
     771    if (largest > 0.001 * fabs(alpha_)) {
    783772      printf("bad\n");
    784773    }
    785774  }
    786775#endif
    787   int numberPivots=abcFactorization_->pivots();
     776  int numberPivots = abcFactorization_->pivots();
    788777  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
    789778  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
     
    791780  if (largestPrimalError_ > 10.0)
    792781    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_))) {
     782  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > checkValue * (1.0 + fabs(alpha_))) {
    795783    handler_->message(CLP_DUAL_CHECK, messages_)
    796784      << btranAlpha_
     
    802790      test = 1.0e-1 * fabs(alpha_);
    803791    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) {
     792      test = 10.0 * checkValue; //1.0e-4 * (1.0 + fabs(alpha_));
     793    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > test);
     794    double derror = CoinMin(fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)), 1.0) * 0.9999e7;
     795    int error = 0;
     796    while (derror > 1.0) {
    810797      error++;
    811798      derror *= 0.1;
    812799    }
    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_)
     800    int frequency[8] = { 99999999, 100, 10, 2, 1, 1, 1, 1 };
     801    int newFactorFrequency = CoinMin(forceFactorization_, frequency[error]);
     802#if ABC_NORMAL_DEBUG > 0
     803    if (newFactorFrequency < forceFactorization_)
    817804      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;
     805        fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)), numberPivots,
     806        forceFactorization_, newFactorFrequency);
     807#endif
     808    if (!numberPivots && fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)) > 0.5e-6
     809      && abcFactorization_->pivotTolerance() < 0.5)
     810      abcFactorization_->saferTolerances(1.0, -1.03);
     811    forceFactorization_ = newFactorFrequency;
    825812
    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_)));
     813#if ABC_NORMAL_DEBUG > 0
     814    if (fabs(btranAlpha_ + alpha_) < 1.0e-5 * (1.0 + fabs(alpha_))) {
     815      printf("diff (%g,%g) %g check %g\n", btranAlpha_, alpha_, fabs(btranAlpha_ - alpha_), checkValue * (1.0 + fabs(alpha_)));
    830816    }
    831817#endif
    832818    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
     819      if (needFlag && numberPivots < 10) {
     820        // need to reject something
     821        assert(sequenceOut_ >= 0);
     822        char x = isColumn(sequenceOut_) ? 'C' : 'R';
     823        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     824          << x << sequenceWithin(sequenceOut_)
     825          << CoinMessageEol;
     826#if ABC_NORMAL_DEBUG > 0
     827        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
     828          btranAlpha_, alpha_, numberPivots);
     829#endif
     830        // Make safer?
     831        abcFactorization_->saferTolerances(1.0, -1.03);
     832        setFlagged(sequenceOut_);
     833        // so can't happen immediately again
     834        sequenceOut_ = -1;
     835        lastBadIteration_ = numberIterations_; // say be more cautious
    850836      }
    851837      // Make safer?
     
    854840      problemStatus_ = -2; // factorize now
    855841      returnCode = -2;
    856       stateOfIteration_=2;
     842      stateOfIteration_ = 2;
    857843    } else {
    858844      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         }
     845        assert(sequenceOut_ >= 0);
     846        // need to reject something
     847        char x = isColumn(sequenceOut_) ? 'C' : 'R';
     848        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     849          << x << sequenceWithin(sequenceOut_)
     850          << CoinMessageEol;
     851#if ABC_NORMAL_DEBUG > 3
     852        printf("flag a %g %g\n", btranAlpha_, alpha_);
     853#endif
     854        setFlagged(sequenceOut_);
     855        // so can't happen immediately again
     856        sequenceOut_ = -1;
     857        //abcProgress_.clearBadTimes();
     858        lastBadIteration_ = numberIterations_; // say be more cautious
     859        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
     860          //printf("I think should declare infeasible\n");
     861          problemStatus_ = 1;
     862          returnCode = 1;
     863          stateOfIteration_ = 2;
     864        } else {
     865          stateOfIteration_ = 1;
     866        }
     867        // Make safer?
     868        if (abcFactorization_->pivotTolerance() < 0.999 && stateOfIteration_ == 1) {
     869          // change tolerance and re-invert
     870          abcFactorization_->saferTolerances(1.0, -1.03);
     871          problemStatus_ = -6; // factorize now
     872          returnCode = -2;
     873          stateOfIteration_ = 2;
     874        }
    889875      }
    890876    }
     
    892878  if (!stateOfIteration_) {
    893879    // check update
    894     //int updateStatus = 
    895 /*
     880    //int updateStatus =
     881    /*
    896882  returns
    897883  0 - OK
     
    901887  5 - take iteration then re-factorize because of memory
    902888 */
    903     int status=checkReplace();
    904     if (status&&!returnCode)
    905       returnCode=status;
     889    int status = checkReplace();
     890    if (status && !returnCode)
     891      returnCode = status;
    906892  }
    907893  //clearArrays(arrayForFlipRhs_);
    908   if (stateOfIteration_&&numberFlipped) {
     894  if (stateOfIteration_ && numberFlipped) {
    909895    //usefulArray_[arrayForTableauRow_].compact();
    910896    // move variables back
     
    915901}
    916902#endif
    917 #if ABC_PARALLEL==2
    918 #if ABC_NORMAL_DEBUG>0
     903#if ABC_PARALLEL == 2
     904#if ABC_NORMAL_DEBUG > 0
    919905// for conflicts
    920 int cilk_conflict=0;
     906int cilk_conflict = 0;
    921907#endif
    922908#ifdef EARLY_FACTORIZE
    923 static 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);
     909static int doEarlyFactorization(AbcSimplexDual *dual)
     910{
     911  int returnCode = cilk_spawn dual->whileIteratingParallel(123456789);
     912  CoinIndexedVector &vector = *dual->usefulArray(ABC_NUMBER_USEFUL - 1);
     913  int status = dual->earlyFactorization()->factorize(dual, vector);
    928914#if 0
    929915  // Is this safe
     
    955941   +3 max iterations
    956942*/
    957 int
    958 AbcSimplexDual::whileIteratingCilk()
     943int AbcSimplexDual::whileIteratingCilk()
    959944{
    960945  /* arrays
     
    976961  int returnCode = -1;
    977962#define DELAYED_UPDATE
    978   arrayForBtran_=0;
     963  arrayForBtran_ = 0;
    979964  setUsedArray(arrayForBtran_);
    980   arrayForFtran_=1;
     965  arrayForFtran_ = 1;
    981966  setUsedArray(arrayForFtran_);
    982   arrayForFlipBounds_=2;
     967  arrayForFlipBounds_ = 2;
    983968  setUsedArray(arrayForFlipBounds_);
    984   arrayForTableauRow_=3;
     969  arrayForTableauRow_ = 3;
    985970  setUsedArray(arrayForTableauRow_);
    986   arrayForDualColumn_=4;
     971  arrayForDualColumn_ = 4;
    987972  setUsedArray(arrayForDualColumn_);
    988   arrayForReplaceColumn_=6; //4;
     973  arrayForReplaceColumn_ = 6; //4;
    989974  setUsedArray(arrayForReplaceColumn_);
    990975#ifndef MOVE_UPDATE_WEIGHTS
    991   arrayForFlipRhs_=5;
     976  arrayForFlipRhs_ = 5;
    992977  setUsedArray(arrayForFlipRhs_);
    993978#else
    994   arrayForFlipRhs_=0;
     979  arrayForFlipRhs_ = 0;
    995980  // use for weights
    996981  setUsedArray(5);
    997982#endif
    998983  dualPivotRow();
    999   lastPivotRow_=pivotRow_;
     984  lastPivotRow_ = pivotRow_;
    1000985  if (pivotRow_ >= 0) {
    1001986    // we found a pivot row
     
    10421027    int numberPivots = abcFactorization_->maximumPivots();
    10431028#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) {
     1029    int numberEarly = 0;
     1030    if (numberPivots > 20 && (numberEarly_ & 0xffff) > 5) {
     1031      numberEarly = numberEarly_ & 0xffff;
     1032      numberPivots = CoinMax(numberPivots - numberEarly - abcFactorization_->pivots(), numberPivots / 2);
     1033    }
     1034    returnCode = whileIteratingParallel(numberPivots);
     1035    if (returnCode == -99) {
    10511036      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
     1037        if (!abcEarlyFactorization_)
     1038          abcEarlyFactorization_ = new AbcSimplexFactorization(*abcFactorization_);
     1039        // create pivot list
     1040        CoinIndexedVector &vector = usefulArray_[ABC_NUMBER_USEFUL - 1];
     1041        vector.checkClear();
     1042        int *indices = vector.getIndices();
     1043        int capacity = vector.capacity();
     1044        int numberNeeded = numberRows_ + 2 * numberEarly + 1;
     1045        if (capacity < numberNeeded) {
     1046          vector.reserve(numberNeeded);
     1047          capacity = vector.capacity();
     1048        }
     1049        int numberBasic = 0;
     1050        for (int i = 0; i < numberRows_; i++) {
     1051          int iSequence = abcPivotVariable_[i];
     1052          if (iSequence < numberRows_)
     1053            indices[numberBasic++] = iSequence;
     1054        }
     1055        vector.setNumElements(numberRows_ - numberBasic);
     1056        for (int i = 0; i < numberRows_; i++) {
     1057          int iSequence = abcPivotVariable_[i];
     1058          if (iSequence >= numberRows_)
     1059            indices[numberBasic++] = iSequence;
     1060        }
     1061        assert(numberBasic == numberRows_);
     1062        indices[capacity - 1] = 0;
     1063        // could set iterations to -1 for safety
    10791064#if 0
    10801065        cilk_spawn doEarlyFactorization(this);
     
    10821067        cilk_sync;
    10831068#else
    1084         returnCode=doEarlyFactorization(this);
     1069        returnCode = doEarlyFactorization(this);
    10851070#endif
    10861071      } else {
    1087         returnCode=whileIteratingParallel(123456789);
     1072        returnCode = whileIteratingParallel(123456789);
    10881073      }
    10891074    }
    10901075#else
    1091     returnCode=whileIteratingParallel(numberPivots);
     1076    returnCode = whileIteratingParallel(numberPivots);
    10921077#endif
    10931078    usefulArray_[arrayForTableauRow_].compact();
     
    10981083    // no pivot row on first try
    10991084    clearArrays(-1);
    1100     returnCode=noPivotRow();
     1085    returnCode = noPivotRow();
    11011086  }
    11021087  //printStuff();
     
    11071092}
    11081093// Create dual pricing vector
    1109 void
    1110 AbcSimplexDual::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;
     1094void AbcSimplexDual::createDualPricingVectorCilk()
     1095{
     1096#if CILK_CONFLICT > 0
     1097  if (cilk_conflict & 15 != 0) {
     1098    printf("cilk_conflict %d!\n", cilk_conflict);
     1099    cilk_conflict = 0;
    11161100  }
    11171101#endif
    11181102#ifndef NDEBUG
    1119 #if ABC_NORMAL_DEBUG>3
    1120     checkArrays();
     1103#if ABC_NORMAL_DEBUG > 3
     1104  checkArrays();
    11211105#endif
    11221106#endif
     
    11381122#if MOVE_REPLACE_PART1A > 0
    11391123  } else {
    1140     cilk_spawn abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
     1124    cilk_spawn abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_], pivotRow_);
    11411125    usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_);
    1142     abcFactorization_->updateColumnTransposeCpu(usefulArray_[arrayForBtran_],1);
     1126    abcFactorization_->updateColumnTransposeCpu(usefulArray_[arrayForBtran_], 1);
    11431127    cilk_sync;
    11441128  }
     
    11461130  sequenceIn_ = -1;
    11471131}
    1148 void
    1149 AbcSimplexDual::getTableauColumnPart1Cilk()
     1132void AbcSimplexDual::getTableauColumnPart1Cilk()
    11501133{
    11511134#if ABC_DEBUG
    11521135  {
    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;
     1136    const double *work = usefulArray_[arrayForTableauRow_].denseVector();
     1137    int number = usefulArray_[arrayForTableauRow_].getNumElements();
     1138    const int *which = usefulArray_[arrayForTableauRow_].getIndices();
     1139    for (int i = 0; i < number; i++) {
     1140      if (which[i] == sequenceIn_) {
     1141        assert(alpha_ == work[i]);
     1142        break;
    11601143      }
    11611144    }
     
    11691152  abcFactorization_->updateColumnFTPart1(usefulArray_[arrayForFtran_]);
    11701153}
    1171 int
    1172 AbcSimplexDual::getTableauColumnFlipAndStartReplaceCilk()
     1154int AbcSimplexDual::getTableauColumnFlipAndStartReplaceCilk()
    11731155{
    11741156  // move checking stuff down into called functions
     
    11821164  cilk_spawn checkReplacePart1();
    11831165#endif
    1184   numberFlipped=flipBounds();
     1166  numberFlipped = flipBounds();
    11851167  cilk_sync;
    11861168#else
    11871169  if (abcFactorization_->usingFT()) {
    11881170    cilk_spawn getTableauColumnPart2();
    1189     ftAlpha_ = cilk_spawn abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
    1190     numberFlipped=flipBounds();
     1171    ftAlpha_ = cilk_spawn abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_], pivotRow_);
     1172    numberFlipped = flipBounds();
    11911173    cilk_sync;
    11921174  } else {
    11931175    cilk_spawn getTableauColumnPart2();
    1194     numberFlipped=flipBounds();
     1176    numberFlipped = flipBounds();
    11951177    cilk_sync;
    11961178  }
     
    12051187    4 - break and say infeasible
    12061188 */
    1207   int returnCode=0;
     1189  int returnCode = 0;
    12081190  // amount primal will move
    12091191  movement_ = -dualOut_ * directionOut_ / alpha_;
    12101192  // see if update stable
    1211 #if ABC_NORMAL_DEBUG>3
     1193#if ABC_NORMAL_DEBUG > 3
    12121194  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_)) {
     1195    double ft = ftAlpha_ * abcFactorization_->pivotRegion()[pivotRow_];
     1196    double diff1 = fabs(alpha_ - btranAlpha_);
     1197    double diff2 = fabs(fabs(alpha_) - fabs(ft));
     1198    double diff3 = fabs(fabs(ft) - fabs(btranAlpha_));
     1199    double largest = CoinMax(CoinMax(diff1, diff2), diff3);
     1200    printf("btran alpha %g, ftran alpha %g ftAlpha %g largest diff %g\n",
     1201      btranAlpha_, alpha_, ft, largest);
     1202    if (largest > 0.001 * fabs(alpha_)) {
    12211203      printf("bad\n");
    12221204    }
    12231205  }
    12241206#endif
    1225   int numberPivots=abcFactorization_->pivots();
     1207  int numberPivots = abcFactorization_->pivots();
    12261208  //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6;
    12271209  double checkValue = numberPivots ? 1.0e-7 : 1.0e-5;
     
    12291211  if (largestPrimalError_ > 10.0)
    12301212    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_))) {
     1213  if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > checkValue * (1.0 + fabs(alpha_))) {
    12331214    handler_->message(CLP_DUAL_CHECK, messages_)
    12341215      << btranAlpha_
     
    12401221      test = 1.0e-1 * fabs(alpha_);
    12411222    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) {
     1223      test = CoinMin(10.0 * checkValue, 1.0e-4 * (1.0 + fabs(alpha_)));
     1224    bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > test);
     1225    double derror = CoinMin(fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)), 1.0) * 0.9999e7;
     1226    int error = 0;
     1227    while (derror > 1.0) {
    12481228      error++;
    12491229      derror *= 0.1;
    12501230    }
    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_)
     1231    int frequency[8] = { 99999999, 100, 10, 2, 1, 1, 1, 1 };
     1232    int newFactorFrequency = CoinMin(forceFactorization_, frequency[error]);
     1233#if ABC_NORMAL_DEBUG > 0
     1234    if (newFactorFrequency < forceFactorization_)
    12551235      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;
     1236        fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)), numberPivots,
     1237        forceFactorization_, newFactorFrequency);
     1238#endif
     1239    if (!numberPivots && fabs(btranAlpha_ - alpha_) / (1.0 + fabs(alpha_)) > 0.5e-6
     1240      && abcFactorization_->pivotTolerance() < 0.5)
     1241      abcFactorization_->saferTolerances(1.0, -1.03);
     1242    forceFactorization_ = newFactorFrequency;
    12631243
    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_)));
     1244#if ABC_NORMAL_DEBUG > 0
     1245    if (fabs(btranAlpha_ + alpha_) < 1.0e-5 * (1.0 + fabs(alpha_))) {
     1246      printf("diff (%g,%g) %g check %g\n", btranAlpha_, alpha_, fabs(btranAlpha_ - alpha_), checkValue * (1.0 + fabs(alpha_)));
    12681247    }
    12691248#endif
    12701249    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
     1250      if (needFlag && numberPivots < 10) {
     1251        // need to reject something
     1252        assert(sequenceOut_ >= 0);
     1253        char x = isColumn(sequenceOut_) ? 'C' : 'R';
     1254        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     1255          << x << sequenceWithin(sequenceOut_)
     1256          << CoinMessageEol;
     1257#if ABC_NORMAL_DEBUG > 0
     1258        printf("flag %d as alpha  %g %g - after %d pivots\n", sequenceOut_,
     1259          btranAlpha_, alpha_, numberPivots);
     1260#endif
     1261        // Make safer?
     1262        abcFactorization_->saferTolerances(1.0, -1.03);
     1263        setFlagged(sequenceOut_);
     1264        // so can't happen immediately again
     1265        sequenceOut_ = -1;
     1266        lastBadIteration_ = numberIterations_; // say be more cautious
    12881267      }
    12891268      // Make safer?
     
    12921271      problemStatus_ = -2; // factorize now
    12931272      returnCode = -2;
    1294       stateOfIteration_=2;
     1273      stateOfIteration_ = 2;
    12951274    } else {
    12961275      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         }
     1276        assert(sequenceOut_ >= 0);
     1277        // need to reject something
     1278        char x = isColumn(sequenceOut_) ? 'C' : 'R';
     1279        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     1280          << x << sequenceWithin(sequenceOut_)
     1281          << CoinMessageEol;
     1282#if ABC_NORMAL_DEBUG > 3
     1283        printf("flag a %g %g\n", btranAlpha_, alpha_);
     1284#endif
     1285        setFlagged(sequenceOut_);
     1286        // so can't happen immediately again
     1287        sequenceOut_ = -1;
     1288        //abcProgress_.clearBadTimes();
     1289        lastBadIteration_ = numberIterations_; // say be more cautious
     1290        if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) {
     1291          //printf("I think should declare infeasible\n");
     1292          problemStatus_ = 1;
     1293          returnCode = 1;
     1294          stateOfIteration_ = 2;
     1295        } else {
     1296          stateOfIteration_ = 1;
     1297        }
     1298        // Make safer?
     1299        if (abcFactorization_->pivotTolerance() < 0.999 && stateOfIteration_ == 1) {
     1300          // change tolerance and re-invert
     1301          abcFactorization_->saferTolerances(1.0, -1.03);
     1302          problemStatus_ = -6; // factorize now
     1303          returnCode = -2;
     1304          stateOfIteration_ = 2;
     1305        }
    13271306      }
    13281307    }
     
    13301309  if (!stateOfIteration_) {
    13311310    // check update
    1332     //int updateStatus = 
    1333 /*
     1311    //int updateStatus =
     1312    /*
    13341313  returns
    13351314  0 - OK
     
    13391318  5 - take iteration then re-factorize because of memory
    13401319 */
    1341     int status=checkReplace();
    1342     if (status&&!returnCode)
    1343       returnCode=status;
     1320    int status = checkReplace();
     1321    if (status && !returnCode)
     1322      returnCode = status;
    13441323  }
    13451324  //clearArrays(arrayForFlipRhs_);
    1346   if (stateOfIteration_&&numberFlipped) {
     1325  if (stateOfIteration_ && numberFlipped) {
    13471326    //usefulArray_[arrayForTableauRow_].compact();
    13481327    // move variables back
     
    13521331  return returnCode;
    13531332}
    1354 int
    1355 AbcSimplexDual::whileIteratingParallel(int numberIterations)
    1356 {
    1357   int returnCode=-1;
     1333int AbcSimplexDual::whileIteratingParallel(int numberIterations)
     1334{
     1335  int returnCode = -1;
    13581336#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;;
     1337  int savePivot = -1;
     1338  CoinIndexedVector &early = usefulArray_[ABC_NUMBER_USEFUL - 1];
     1339  int *pivotIndices = early.getIndices();
     1340  double *pivotValues = early.denseVector();
     1341  int pivotCountPosition = early.capacity() - 1;
     1342  int numberSaved = 0;
     1343  if (numberIterations == 123456789)
     1344    savePivot = pivotCountPosition;
     1345  ;
    13671346#endif
    13681347  numberIterations += numberIterations_;
     
    13771356      2 - don't take and break
    13781357    */
    1379     stateOfIteration_=0;
    1380     returnCode=-1;
     1358    stateOfIteration_ = 0;
     1359    returnCode = -1;
    13811360    // put row of tableau in usefulArray[arrayForTableauRow_]
    13821361    /*
     
    13871366    */
    13881367    //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) {
     1368    double saveAcceptable = acceptablePivot_;
     1369    if (largestPrimalError_ > 1.0e-5 || largestDualError_ > 1.0e-5) {
     1370      //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) {
    13921371      if (!abcFactorization_->pivots())
    1393         acceptablePivot_*=1.0e2;
    1394       else if (abcFactorization_->pivots()<5)
    1395         acceptablePivot_*=1.0e1;
     1372        acceptablePivot_ *= 1.0e2;
     1373      else if (abcFactorization_->pivots() < 5)
     1374        acceptablePivot_ *= 1.0e1;
    13961375    }
    13971376#ifdef MOVE_UPDATE_WEIGHTS
    1398     // copy btran across 
     1377    // copy btran across
    13991378    usefulArray_[5].copy(usefulArray_[arrayForBtran_]);
    1400     cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[5]);;
     1379    cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[5]);
     1380    ;
    14011381#endif
    14021382    dualColumn1();
    1403     acceptablePivot_=saveAcceptable;
    1404     if ((stateOfProblem_&VALUES_PASS)!=0) {
     1383    acceptablePivot_ = saveAcceptable;
     1384    if ((stateOfProblem_ & VALUES_PASS) != 0) {
    14051385      // see if can just move dual
    1406       if (fabs(upperTheta_-fabs(abcDj_[sequenceOut_]))<dualTolerance_) {
    1407         stateOfIteration_=-1;
     1386      if (fabs(upperTheta_ - fabs(abcDj_[sequenceOut_])) < dualTolerance_) {
     1387        stateOfIteration_ = -1;
    14081388      }
    14091389    }
    14101390    if (!stateOfIteration_) {
    14111391#ifndef MOVE_UPDATE_WEIGHTS
    1412       cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[arrayForBtran_]);;
     1392      cilk_spawn abcDualRowPivot_->updateWeightsOnly(usefulArray_[arrayForBtran_]);
     1393      ;
    14131394#endif
    14141395      // get sequenceIn_
     
    14161397      //stateOfIteration_=0;
    14171398      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;
     1399        // normal iteration
     1400        // update the incoming column
     1401        //arrayForReplaceColumn_=getAvailableArray();
     1402        getTableauColumnFlipAndStartReplaceCilk();
     1403        //usleep(1000);
     1404      } else if (stateOfIteration_ != -1) {
     1405        stateOfIteration_ = 2;
    14251406      }
    14261407    }
     
    14301411      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
    14311412      if (status >= 0) {
    1432         problemStatus_ = 5;
    1433         secondaryStatus_ = ClpEventHandler::endOfIteration;
    1434         returnCode = 4;
    1435         stateOfIteration_=2;
     1413        problemStatus_ = 5;
     1414        secondaryStatus_ = ClpEventHandler::endOfIteration;
     1415        returnCode = 4;
     1416        stateOfIteration_ = 2;
    14361417      }
    14371418    }
    14381419    if (!stateOfIteration_) {
    1439       //        assert (stateOfIteration_!=1); 
    1440       int whatNext=0;
     1420      //        assert (stateOfIteration_!=1);
     1421      int whatNext = 0;
    14411422      whatNext = housekeeping();
    14421423#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++;
     1424      if (savePivot >= 0) {
     1425        // save pivot
     1426        pivotIndices[--savePivot] = sequenceOut_;
     1427        pivotValues[savePivot] = alpha_;
     1428        pivotIndices[--savePivot] = sequenceIn_;
     1429        pivotValues[savePivot] = btranAlpha_;
     1430        numberSaved++;
    14501431      }
    14511432#endif
    14521433      if (whatNext == 1) {
    1453         problemStatus_ = -2; // refactorize
    1454         usefulArray_[arrayForTableauRow_].compact();
     1434        problemStatus_ = -2; // refactorize
     1435        usefulArray_[arrayForTableauRow_].compact();
    14551436      } else if (whatNext == 2) {
    1456         // maximum iterations or equivalent
    1457         problemStatus_ = 3;
    1458         returnCode = 3;
    1459         stateOfIteration_=2;
     1437        // maximum iterations or equivalent
     1438        problemStatus_ = 3;
     1439        returnCode = 3;
     1440        stateOfIteration_ = 2;
    14601441      }
    14611442#ifdef EARLY_FACTORIZE
    1462       if (savePivot>=0) {
    1463         // tell worker can update this
    1464         pivotIndices[pivotCountPosition]=numberSaved;
     1443      if (savePivot >= 0) {
     1444        // tell worker can update this
     1445        pivotIndices[pivotCountPosition] = numberSaved;
    14651446      }
    14661447#endif
    14671448    } else {
    14681449      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;
     1450      if (stateOfIteration_ < 0) {
     1451        // move dual in dual values pass
     1452        theta_ = abcDj_[sequenceOut_];
     1453        updateDualsInDual();
     1454        abcDj_[sequenceOut_] = 0.0;
     1455        sequenceOut_ = -1;
    14751456      }
    14761457      // clear all arrays
     
    14781459    }
    14791460    // at this stage sequenceIn_ may be <0
    1480     if (sequenceIn_<0&&sequenceOut_>=0) {
     1461    if (sequenceIn_ < 0 && sequenceOut_ >= 0) {
    14811462      usefulArray_[arrayForTableauRow_].compact();
    14821463      // no incoming column is valid
    1483       returnCode=noPivotColumn();
    1484     }
    1485     if (stateOfIteration_==2) {
     1464      returnCode = noPivotColumn();
     1465    }
     1466    if (stateOfIteration_ == 2) {
    14861467      usefulArray_[arrayForTableauRow_].compact();
    1487       sequenceIn_=-1;
     1468      sequenceIn_ = -1;
    14881469#ifdef ABC_LONG_FACTORIZATION
    14891470      abcFactorization_->clearHiddenArrays();
     
    14911472      break;
    14921473    }
    1493     if (problemStatus_!=-1) {
     1474    if (problemStatus_ != -1) {
    14941475#ifndef MOVE_UPDATE_WEIGHTS
    1495       abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_]);
     1476      abcDualRowPivot_->updateWeights2(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_]);
    14961477#else
    1497       abcDualRowPivot_->updateWeights2(usefulArray_[5],usefulArray_[arrayForFtran_]);
     1478      abcDualRowPivot_->updateWeights2(usefulArray_[5], usefulArray_[arrayForFtran_]);
    14981479#endif
    14991480#ifdef ABC_LONG_FACTORIZATION
     
    15031484      break;
    15041485    }
    1505     if (stateOfIteration_==1) {
     1486    if (stateOfIteration_ == 1) {
    15061487      // clear all arrays
    15071488      clearArrays(-2);
     
    15101491#endif
    15111492    }
    1512     if (stateOfIteration_==0) {
     1493    if (stateOfIteration_ == 0) {
    15131494      // can do these in parallel
    15141495      // No idea why I need this - but otherwise runs not repeatable (try again??)
     
    15171498      int lastSequenceOut;
    15181499      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;
     1500      if (firstFree_ < 0) {
     1501        // can do in parallel
     1502        cilk_spawn replaceColumnPart3();
     1503        updatePrimalSolution();
     1504        swapPrimalStuff();
     1505        // dualRow will go to virtual row pivot choice algorithm
     1506        // make sure values pass off if it should be
     1507        // use Btran array and clear inside dualPivotRow (if used)
     1508        lastSequenceOut = sequenceOut_;
     1509        lastDirectionOut = directionOut_;
     1510        dualPivotRow();
     1511        cilk_sync;
    15311512      } 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_;
     1513        // be more careful as dualPivotRow may do update
     1514        cilk_spawn replaceColumnPart3();
     1515        updatePrimalSolution();
     1516        swapPrimalStuff();
     1517        // dualRow will go to virtual row pivot choice algorithm
     1518        // make sure values pass off if it should be
     1519        // use Btran array and clear inside dualPivotRow (if used)
     1520        lastSequenceOut = sequenceOut_;
     1521        lastDirectionOut = directionOut_;
     1522        cilk_sync;
     1523        dualPivotRow();
     1524      }
     1525      lastPivotRow_ = pivotRow_;
    15451526      if (pivotRow_ >= 0) {
    1546         // we found a pivot row
    1547         // create as packed
    1548         // dual->checkReplacePart1a();
    1549         createDualPricingVectorCilk();
    1550         swapDualStuff(lastSequenceOut,lastDirectionOut);
     1527        // we found a pivot row
     1528        // create as packed
     1529        // dual->checkReplacePart1a();
     1530        createDualPricingVectorCilk();
     1531        swapDualStuff(lastSequenceOut, lastDirectionOut);
    15511532      }
    15521533      cilk_sync;
     
    15541535      // after moving dual in values pass
    15551536      dualPivotRow();
    1556       lastPivotRow_=pivotRow_;
     1537      lastPivotRow_ = pivotRow_;
    15571538      if (pivotRow_ >= 0) {
    1558         // we found a pivot row
    1559         // create as packed
    1560         createDualPricingVectorCilk();
    1561       }
    1562     }
    1563     if (pivotRow_<0) {
     1539        // we found a pivot row
     1540        // create as packed
     1541        createDualPricingVectorCilk();
     1542      }
     1543    }
     1544    if (pivotRow_ < 0) {
    15641545      // no pivot row
    15651546      clearArrays(-1);
    1566       returnCode=noPivotRow();
     1547      returnCode = noPivotRow();
    15671548      break;
    15681549    }
    1569     if (numberIterations == numberIterations_&&problemStatus_==-1) {
    1570       returnCode=-99;
     1550    if (numberIterations == numberIterations_ && problemStatus_ == -1) {
     1551      returnCode = -99;
    15711552      break;
    15721553    }
     
    15931574#endif
    15941575#endif
    1595 void
    1596 AbcSimplexDual::updatePrimalSolution()
     1576void AbcSimplexDual::updatePrimalSolution()
    15971577{
    15981578  // finish doing weights
    15991579#ifndef MOVE_UPDATE_WEIGHTS
    1600   abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_],usefulArray_[arrayForFtran_],
    1601                                         movement_);
     1580  abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[arrayForBtran_], usefulArray_[arrayForFtran_],
     1581    movement_);
    16021582#else
    1603   abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[5],usefulArray_[arrayForFtran_],
    1604                                         movement_);
     1583  abcDualRowPivot_->updatePrimalSolutionAndWeights(usefulArray_[5], usefulArray_[arrayForFtran_],
     1584    movement_);
    16051585#endif
    16061586}
    16071587int xxInfo[6][8];
    1608 double parallelDual4(AbcSimplexDual * dual)
    1609 {
    1610   int maximumRows=dual->maximumAbcNumberRows();
     1588double parallelDual4(AbcSimplexDual *dual)
     1589{
     1590  int maximumRows = dual->maximumAbcNumberRows();
    16111591  //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();
     1592  CoinIndexedVector &update = *dual->usefulArray(dual->arrayForBtran());
     1593  CoinPartitionedVector &tableauRow = *dual->usefulArray(dual->arrayForTableauRow());
     1594  CoinPartitionedVector &candidateList = *dual->usefulArray(dual->arrayForDualColumn());
     1595  AbcMatrix *matrix = dual->abcMatrix();
    16161596  // do slacks first (move before sync)
    1617   double upperTheta=dual->upperTheta();
     1597  double upperTheta = dual->upperTheta();
    16181598  //assert (upperTheta>-100*dual->dualTolerance()||dual->sequenceIn()>=0||dual->lastFirstFree()>=0);
    16191599#if ABC_PARALLEL
    16201600#define NUMBER_BLOCKS NUMBER_ROW_BLOCKS
    1621   int numberBlocks=CoinMin(NUMBER_BLOCKS,dual->numberCpus());
     1601  int numberBlocks = CoinMin(NUMBER_BLOCKS, dual->numberCpus());
    16221602#else
    16231603#define NUMBER_BLOCKS 1
    1624   int numberBlocks=NUMBER_BLOCKS;
     1604  int numberBlocks = NUMBER_BLOCKS;
    16251605#endif
    16261606#if ABC_PARALLEL
    1627   if (dual->parallelMode()==0)
    1628     numberBlocks=1;
     1607  if (dual->parallelMode() == 0)
     1608    numberBlocks = 1;
    16291609#endif
    16301610  // see if worth using row copy
    16311611  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
     1612  int number = update.getNumElements();
     1613  assert(number);
     1614  bool useRowCopy = (gotRowCopy && (number << 2) < maximumRows);
     1615  assert(numberBlocks == matrix->numberRowBlocks());
     1616#if ABC_PARALLEL == 1
    16371617  // redo so can pass info in with void *
    1638   CoinThreadInfo * infoP = dual->threadInfoPointer();
    1639   int cpuMask=((1<<dual->parallelMode())-1);
    1640   cpuMask += cpuMask<<5;
     1618  CoinThreadInfo *infoP = dual->threadInfoPointer();
     1619  int cpuMask = ((1 << dual->parallelMode()) - 1);
     1620  cpuMask += cpuMask << 5;
    16411621  dual->setStopStart(cpuMask);
    16421622#endif
    16431623  CoinThreadInfo info[NUMBER_BLOCKS];
    1644   const int * starts;
     1624  const int *starts;
    16451625  if (useRowCopy)
    1646     starts=matrix->blockStart();
     1626    starts = matrix->blockStart();
    16471627  else
    1648     starts=matrix->startColumnBlock();
    1649   tableauRow.setPartitions(numberBlocks,starts);
    1650   candidateList.setPartitions(numberBlocks,starts);
     1628    starts = matrix->startColumnBlock();
     1629  tableauRow.setPartitions(numberBlocks, starts);
     1630  candidateList.setPartitions(numberBlocks, starts);
    16511631  //printf("free sequence in %d\n",dual->freeSequenceIn());
    16521632  if (useRowCopy) {
    16531633    // using row copy
    16541634#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;
     1635    if (numberBlocks > 1) {
     1636#if ABC_PARALLEL == 2
     1637      for (int i = 0; i < numberBlocks; i++) {
     1638        info[i].stuff[1] = i;
     1639        info[i].stuff[2] = -1;
     1640        info[i].result = upperTheta;
     1641        info[i].result = cilk_spawn
     1642                           matrix->dualColumn1Row(info[i].stuff[1], COIN_DBL_MAX, info[i].stuff[2],
     1643                             update, tableauRow, candidateList);
     1644      }
     1645      cilk_sync;
    16661646#else
    16671647      // 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];
     1648      for (int i = 0; i < numberBlocks; i++) {
     1649        info[i].status = 5;
     1650        info[i].stuff[1] = i;
     1651        info[i].result = upperTheta;
     1652        if (i < numberBlocks - 1) {
     1653          infoP[i] = info[i];
     1654          if (i == numberBlocks - 2)
     1655            dual->startParallelStuff(5);
     1656        }
     1657      }
     1658      info[numberBlocks - 1].stuff[1] = numberBlocks - 1;
     1659      info[numberBlocks - 1].stuff[2] = -1;
     1660      //double freeAlpha;
     1661      info[numberBlocks - 1].result = matrix->dualColumn1Row(info[numberBlocks - 1].stuff[1],
     1662        upperTheta, info[numberBlocks - 1].stuff[2],
     1663        update, tableauRow, candidateList);
     1664      dual->stopParallelStuff(5);
     1665      for (int i = 0; i < numberBlocks - 1; i++)
     1666        info[i] = infoP[i];
    16881667#endif
    16891668    } else {
    16901669#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;
     1670      info[0].status = 5;
     1671      info[0].stuff[1] = 0;
     1672      info[0].result = upperTheta;
     1673      info[0].stuff[1] = 0;
     1674      info[0].stuff[2] = -1;
    16961675      // worth using row copy
    16971676      //assert (number>2);
    1698       info[0].result =
    1699       matrix->dualColumn1Row(info[0].stuff[1],
    1700                              upperTheta,info[0].stuff[2],
    1701                              update,tableauRow,candidateList);
     1677      info[0].result = matrix->dualColumn1Row(info[0].stuff[1],
     1678        upperTheta, info[0].stuff[2],
     1679        update, tableauRow, candidateList);
    17021680#if ABC_PARALLEL
    17031681    }
     
    17051683  } else { // end of useRowCopy
    17061684#if ABC_PARALLEL
    1707     if (numberBlocks>1) {
    1708 #if ABC_PARALLEL==2
     1685    if (numberBlocks > 1) {
     1686#if ABC_PARALLEL == 2
    17091687      // 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);
     1688      for (int i = 0; i < numberBlocks; i++) {
     1689        info[i].stuff[1] = i;
     1690        info[i].result = upperTheta;
     1691        cilk_spawn
     1692          matrix->dualColumn1Part(info[i].stuff[1], info[i].stuff[2],
     1693            info[i].result,
     1694            update, tableauRow, candidateList);
    17171695      }
    17181696      cilk_sync;
     
    17201698      // parallel 1
    17211699      // 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);
     1700      for (int i = 0; i < numberBlocks; i++) {
     1701        info[i].status = 4;
     1702        info[i].stuff[1] = i;
     1703        info[i].result = upperTheta;
     1704        if (i < numberBlocks - 1) {
     1705          infoP[i] = info[i];
     1706          if (i == numberBlocks - 2)
     1707            dual->startParallelStuff(4);
     1708        }
     1709      }
     1710      matrix->dualColumn1Part(info[numberBlocks - 1].stuff[1], info[numberBlocks - 1].stuff[2],
     1711        info[numberBlocks - 1].result,
     1712        update, tableauRow, candidateList);
    17351713      dual->stopParallelStuff(4);
    1736       for (int i=0;i<numberBlocks-1;i++)
    1737         info[i]=infoP[i];
     1714      for (int i = 0; i < numberBlocks - 1; i++)
     1715        info[i] = infoP[i];
    17381716#endif
    17391717    } else {
    17401718#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);
     1719      info[0].status = 4;
     1720      info[0].stuff[1] = 0;
     1721      info[0].result = upperTheta;
     1722      info[0].stuff[1] = 0;
     1723      matrix->dualColumn1Part(info[0].stuff[1], info[0].stuff[2],
     1724        info[0].result,
     1725        update, tableauRow, candidateList);
    17481726#if ABC_PARALLEL
    17491727    }
     
    17511729  }
    17521730  int sequenceIn[NUMBER_BLOCKS];
    1753   bool anyFree=false;
     1731  bool anyFree = false;
    17541732#if 0
    17551733  if (useRowCopy) {
     
    17751753  }
    17761754#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);
     1755  for (int i = 0; i < numberBlocks; i++) {
     1756    sequenceIn[i] = info[i].stuff[2];
     1757    if (sequenceIn[i] >= 0)
     1758      anyFree = true;
     1759    upperTheta = CoinMin(info[i].result, upperTheta);
    17821760    //assert (info[i].result>-100*dual->dualTolerance()||sequenceIn[i]>=0||dual->lastFirstFree()>=0);
    17831761  }
    17841762  if (anyFree) {
    1785     int *  COIN_RESTRICT index = tableauRow.getIndices();
    1786     double *  COIN_RESTRICT array = tableauRow.denseVector();
     1763    int *COIN_RESTRICT index = tableauRow.getIndices();
     1764    double *COIN_RESTRICT array = tableauRow.denseVector();
    17871765    // 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) {
     1766    double bestFree = 0.0;
     1767    int bestSequence = dual->sequenceIn();
     1768    if (bestSequence >= 0)
     1769      bestFree = dual->alpha();
     1770    for (int i = 0; i < numberBlocks; i++) {
     1771      int iLook = sequenceIn[i];
     1772      if (iLook >= 0) {
     1773        // free variable - search
     1774        int start = starts[i];
     1775        int end = start + tableauRow.getNumElements(i);
     1776        for (int k = start; k < end; k++) {
     1777          if (iLook == index[k]) {
     1778            if (fabs(bestFree) < fabs(array[k])) {
     1779              bestFree = array[k];
     1780              bestSequence = iLook;
     1781            }
     1782            break;
     1783          }
     1784        }
     1785      }
     1786    }
     1787    if (bestSequence >= 0) {
    18101788      double oldValue = dual->djRegion()[bestSequence];
    18111789      dual->setSequenceIn(bestSequence);
     
    18161794  //tableauRow.compact();
    18171795  //candidateList.compact();
    1818 #if 0//ndef NDEBUG
     1796#if 0 //ndef NDEBUG
    18191797  tableauRow.setPackedMode(true);
    18201798  tableauRow.sortPacked();
     
    18241802  return upperTheta;
    18251803}
    1826 #if ABC_PARALLEL==2
     1804#if ABC_PARALLEL == 2
    18271805static void
    1828 parallelDual5a(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) {
     1806parallelDual5a(AbcSimplexFactorization *factorization,
     1807  CoinIndexedVector *whichVector,
     1808  int numberCpu,
     1809  int whichCpu,
     1810  double *weights)
     1811{
     1812  double *COIN_RESTRICT array = whichVector->denseVector();
     1813  int *COIN_RESTRICT which = whichVector->getIndices();
     1814  int numberRows = factorization->numberRows();
     1815  for (int i = whichCpu; i < numberRows; i += numberCpu) {
    18381816    double value = 0.0;
    18391817    array[i] = 1.0;
    18401818    which[0] = i;
    18411819    whichVector->setNumElements(1);
    1842     factorization->updateColumnTransposeCpu(*whichVector,whichCpu);
     1820    factorization->updateColumnTransposeCpu(*whichVector, whichCpu);
    18431821    int number = whichVector->getNumElements();
    18441822    for (int j = 0; j < number; j++) {
    1845       int k= which[j];
     1823      int k = which[j];
    18461824      value += array[k] * array[k];
    18471825      array[k] = 0.0;
     
    18521830}
    18531831#endif
    1854 #if ABC_PARALLEL==2
    1855 void
    1856 parallelDual5(AbcSimplexFactorization * factorization,
    1857               CoinIndexedVector ** whichVector,
    1858               int numberCpu,
    1859               int whichCpu,
    1860               double * weights)
     1832#if ABC_PARALLEL == 2
     1833void parallelDual5(AbcSimplexFactorization *factorization,
     1834  CoinIndexedVector **whichVector,
     1835  int numberCpu,
     1836  int whichCpu,
     1837  double *weights)
    18611838{
    18621839  if (whichCpu) {
    1863     cilk_spawn parallelDual5(factorization,whichVector,numberCpu,whichCpu-1,weights);
    1864     parallelDual5a(factorization,whichVector[whichCpu],numberCpu,whichCpu,weights);
     1840    cilk_spawn parallelDual5(factorization, whichVector, numberCpu, whichCpu - 1, weights);
     1841    parallelDual5a(factorization, whichVector[whichCpu], numberCpu, whichCpu, weights);
    18651842    cilk_sync;
    18661843  } else {
    1867     parallelDual5a(factorization,whichVector[whichCpu],numberCpu,whichCpu,weights);
     1844    parallelDual5a(factorization, whichVector[whichCpu], numberCpu, whichCpu, weights);
    18681845  }
    18691846}
     
    18711848// cilk seems a bit fragile
    18721849#define CILK_FRAGILE 1
    1873 #if CILK_FRAGILE>1
     1850#if CILK_FRAGILE > 1
    18741851#undef cilk_spawn
    18751852#undef cilk_sync
     
    18771854#define cilk_sync
    18781855#define ONWARD 0
    1879 #elif CILK_FRAGILE==1
     1856#elif CILK_FRAGILE == 1
    18801857#define ONWARD 0
    18811858#else
     
    18921869   3 Left Upper Transpose NonUnit
    18931870*/
    1894 static 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);
     1871static void CoinAbcDtrsmFactor(int m, int n, double *COIN_RESTRICT a, int lda)
     1872{
     1873  assert((m & (BLOCKING8 - 1)) == 0 && (n & (BLOCKING8 - 1)) == 0);
     1874  assert(m == BLOCKING8);
    18981875  // 0 Left Lower NoTranspose Unit
    18991876  /* entry for column j and row i (when multiple of BLOCKING8)
    19001877     is at aBlocked+j*m+i*BLOCKING8
    19011878  */
    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++) {
     1879  double *COIN_RESTRICT aBase2 = a;
     1880  double *COIN_RESTRICT bBase2 = aBase2 + lda * BLOCKING8;
     1881  for (int jj = 0; jj < n; jj += BLOCKING8) {
     1882    double *COIN_RESTRICT bBase = bBase2;
     1883    for (int j = jj; j < jj + BLOCKING8; j++) {
    19071884#if 0
    19081885      double * COIN_RESTRICT aBase = aBase2;
     
    19181895#else
    19191896      // unroll - stay in registers - don't test for zeros
    1920       assert (BLOCKING8==8);
     1897      assert(BLOCKING8 == 8);
    19211898      double bValue0 = bBase[0];
    19221899      double bValue1 = bBase[1];
     
    19271904      double bValue6 = bBase[6];
    19281905      double bValue7 = bBase[7];
    1929       bValue1-=bValue0*a[1+0*BLOCKING8];
     1906      bValue1 -= bValue0 * a[1 + 0 * BLOCKING8];
    19301907      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];
     1908      bValue2 -= bValue0 * a[2 + 0 * BLOCKING8];
     1909      bValue3 -= bValue0 * a[3 + 0 * BLOCKING8];
     1910      bValue4 -= bValue0 * a[4 + 0 * BLOCKING8];
     1911      bValue5 -= bValue0 * a[5 + 0 * BLOCKING8];
     1912      bValue6 -= bValue0 * a[6 + 0 * BLOCKING8];
     1913      bValue7 -= bValue0 * a[7 + 0 * BLOCKING8];
     1914      bValue2 -= bValue1 * a[2 + 1 * BLOCKING8];
    19381915      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];
     1916      bValue3 -= bValue1 * a[3 + 1 * BLOCKING8];
     1917      bValue4 -= bValue1 * a[4 + 1 * BLOCKING8];
     1918      bValue5 -= bValue1 * a[5 + 1 * BLOCKING8];
     1919      bValue6 -= bValue1 * a[6 + 1 * BLOCKING8];
     1920      bValue7 -= bValue1 * a[7 + 1 * BLOCKING8];
     1921      bValue3 -= bValue2 * a[3 + 2 * BLOCKING8];
    19451922      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];
     1923      bValue4 -= bValue2 * a[4 + 2 * BLOCKING8];
     1924      bValue5 -= bValue2 * a[5 + 2 * BLOCKING8];
     1925      bValue6 -= bValue2 * a[6 + 2 * BLOCKING8];
     1926      bValue7 -= bValue2 * a[7 + 2 * BLOCKING8];
     1927      bValue4 -= bValue3 * a[4 + 3 * BLOCKING8];
    19511928      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];
     1929      bValue5 -= bValue3 * a[5 + 3 * BLOCKING8];
     1930      bValue6 -= bValue3 * a[6 + 3 * BLOCKING8];
     1931      bValue7 -= bValue3 * a[7 + 3 * BLOCKING8];
     1932      bValue5 -= bValue4 * a[5 + 4 * BLOCKING8];
    19561933      bBase[5] = bValue5;
    1957       bValue6-=bValue4*a[6+4*BLOCKING8];
    1958       bValue7-=bValue4*a[7+4*BLOCKING8];
    1959       bValue6-=bValue5*a[6+5*BLOCKING8];
     1934      bValue6 -= bValue4 * a[6 + 4 * BLOCKING8];
     1935      bValue7 -= bValue4 * a[7 + 4 * BLOCKING8];
     1936      bValue6 -= bValue5 * a[6 + 5 * BLOCKING8];
    19601937      bBase[6] = bValue6;
    1961       bValue7-=bValue5*a[7+5*BLOCKING8];
    1962       bValue7-=bValue6*a[7+6*BLOCKING8];
     1938      bValue7 -= bValue5 * a[7 + 5 * BLOCKING8];
     1939      bValue7 -= bValue6 * a[7 + 6 * BLOCKING8];
    19631940      bBase[7] = bValue7;
    19641941#endif
    19651942      bBase += BLOCKING8;
    19661943    }
    1967     bBase2 +=lda*BLOCKING8;
     1944    bBase2 += lda * BLOCKING8;
    19681945  }
    19691946}
     
    19711948#define CILK_DTRSM 32
    19721949static 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);
     1950  int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     1951{
     1952  int mm = CoinMin(kkk + UNROLL_DTRSM * BLOCKING8, m);
     1953  assert((last - first) % BLOCKING8 == 0);
     1954  if (last - first > CILK_DTRSM) {
     1955    int mid = ((first + last) >> 4) << 3;
     1956    cilk_spawn dtrsm0(kkk, first, mid, m, a, b);
     1957    dtrsm0(kkk, mid, last, m, a, b);
    19811958    cilk_sync;
    19821959  } 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) {
     1960    const double *COIN_RESTRICT aBaseA = a + UNROLL_DTRSM * BLOCKING8X8 + kkk * BLOCKING8;
     1961    aBaseA += (first - mm) * BLOCKING8 - BLOCKING8X8;
     1962    aBaseA += m * kkk;
     1963    for (int ii = first; ii < last; ii += BLOCKING8) {
    19871964      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
     1965      const double *COIN_RESTRICT aBaseB = aBaseA;
     1966      double *COIN_RESTRICT bBaseI = b + ii;
     1967      for (int kk = kkk; kk < mm; kk += BLOCKING8) {
     1968        double *COIN_RESTRICT bBase = b + kk;
     1969        const double *COIN_RESTRICT aBase2 = aBaseB; //a+UNROLL_DTRSM*BLOCKING8X8+m*kk+kkk*BLOCKING8;
     1970        //aBase2 += (ii-mm)*BLOCKING8;
     1971        //assert (aBase2==aBaseB);
     1972        aBaseB += m * BLOCKING8;
     1973#if AVX2 != 2
    19971974#define ALTERNATE_INNER
    19981975#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         }
     1976        for (int k = 0; k < BLOCKING8; k++) {
     1977          //coin_prefetch_const(aBase2+BLOCKING8);
     1978          for (int i = 0; i < BLOCKING8; i++) {
     1979            bBaseI[i] -= bBase[k] * aBase2[i];
     1980          }
     1981          aBase2 += BLOCKING8;
     1982        }
    20061983#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;
     1984        double b0 = bBaseI[0];
     1985        double b1 = bBaseI[1];
     1986        double b2 = bBaseI[2];
     1987        double b3 = bBaseI[3];
     1988        double b4 = bBaseI[4];
     1989        double b5 = bBaseI[5];
     1990        double b6 = bBaseI[6];
     1991        double b7 = bBaseI[7];
     1992        for (int k = 0; k < BLOCKING8; k++) {
     1993          //coin_prefetch_const(aBase2+BLOCKING8);
     1994          double bValue = bBase[k];
     1995          b0 -= bValue * aBase2[0];
     1996          b1 -= bValue * aBase2[1];
     1997          b2 -= bValue * aBase2[2];
     1998          b3 -= bValue * aBase2[3];
     1999          b4 -= bValue * aBase2[4];
     2000          b5 -= bValue * aBase2[5];
     2001          b6 -= bValue * aBase2[6];
     2002          b7 -= bValue * aBase2[7];
     2003          aBase2 += BLOCKING8;
     2004        }
     2005        bBaseI[0] = b0;
     2006        bBaseI[1] = b1;
     2007        bBaseI[2] = b2;
     2008        bBaseI[3] = b3;
     2009        bBaseI[4] = b4;
     2010        bBaseI[5] = b5;
     2011        bBaseI[6] = b6;
     2012        bBaseI[7] = b7;
    20362013#endif
    20372014#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 }
    2058 void CoinAbcDtrsm0(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
    2059 {
    2060   assert ((m&(BLOCKING8-1))==0);
     2015        __m256d b0 = _mm256_load_pd(bBaseI);
     2016        __m256d b1 = _mm256_load_pd(bBaseI + 4);
     2017        for (int k = 0; k < BLOCKING8; k++) {
     2018          //coin_prefetch_const(aBase2+BLOCKING8);
     2019          __m256d bb = _mm256_broadcast_sd(bBase + k);
     2020          __m256d a0 = _mm256_load_pd(aBase2);
     2021          b0 -= a0 * bb;
     2022          __m256d a1 = _mm256_load_pd(aBase2 + 4);
     2023          b1 -= a1 * bb;
     2024          aBase2 += BLOCKING8;
     2025        }
     2026        //_mm256_store_pd ((bBaseI), (b0));
     2027        *reinterpret_cast< __m256d * >(bBaseI) = b0;
     2028        //_mm256_store_pd (bBaseI+4, b1);
     2029        *reinterpret_cast< __m256d * >(bBaseI + 4) = b1;
     2030#endif
     2031      }
     2032    }
     2033  }
     2034}
     2035void CoinAbcDtrsm0(int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     2036{
     2037  assert((m & (BLOCKING8 - 1)) == 0);
    20612038  // 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);
     2039  for (int kkk = 0; kkk < m; kkk += UNROLL_DTRSM * BLOCKING8) {
     2040    int mm = CoinMin(m, kkk + UNROLL_DTRSM * BLOCKING8);
     2041    for (int kk = kkk; kk < mm; kk += BLOCKING8) {
     2042      const double *COIN_RESTRICT aBase2 = a + kk * (m + BLOCKING8);
     2043      double *COIN_RESTRICT bBase = b + kk;
     2044      for (int k = 0; k < BLOCKING8; k++) {
     2045        for (int i = k + 1; i < BLOCKING8; i++) {
     2046          bBase[i] -= bBase[k] * aBase2[i];
     2047        }
     2048        aBase2 += BLOCKING8;
     2049      }
     2050      for (int ii = kk + BLOCKING8; ii < mm; ii += BLOCKING8) {
     2051        double *COIN_RESTRICT bBaseI = b + ii;
     2052        for (int k = 0; k < BLOCKING8; k++) {
     2053          //coin_prefetch_const(aBase2+BLOCKING8);
     2054          for (int i = 0; i < BLOCKING8; i++) {
     2055            bBaseI[i] -= bBase[k] * aBase2[i];
     2056          }
     2057          aBase2 += BLOCKING8;
     2058        }
     2059      }
     2060    }
     2061    dtrsm0(kkk, mm, m, m, a, b);
    20852062  }
    20862063}
    20872064static 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);
     2065  int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     2066{
     2067  int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2068  assert((last - first) % BLOCKING8 == 0);
     2069  if (last - first > CILK_DTRSM) {
     2070    int mid = ((first + last) >> 4) << 3;
     2071    cilk_spawn dtrsm1(kkk, first, mid, m, a, b);
     2072    dtrsm1(kkk, mid, last, m, a, b);
    20962073    cilk_sync;
    20972074  } 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
     2075    for (int iii = last - BLOCKING8; iii >= first; iii -= BLOCKING8) {
     2076      double *COIN_RESTRICT bBase2 = b + iii;
     2077      const double *COIN_RESTRICT aBaseA = a + BLOCKING8X8 + BLOCKING8 * iii;
     2078      for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2079        double *COIN_RESTRICT bBase = b + ii;
     2080        const double *COIN_RESTRICT aBase = aBaseA + m * ii;
     2081#if AVX2 != 2
    21062082#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         }
     2083        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2084          aBase -= BLOCKING8;
     2085          //coin_prefetch_const(aBase-BLOCKING8);
     2086          for (int k = BLOCKING8 - 1; k >= 0; k--) {
     2087            bBase2[k] -= bBase[i] * aBase[k];
     2088          }
     2089        }
    21142090#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;
     2091        double b0 = bBase2[0];
     2092        double b1 = bBase2[1];
     2093        double b2 = bBase2[2];
     2094        double b3 = bBase2[3];
     2095        double b4 = bBase2[4];
     2096        double b5 = bBase2[5];
     2097        double b6 = bBase2[6];
     2098        double b7 = bBase2[7];
     2099        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2100          aBase -= BLOCKING8;
     2101          //coin_prefetch_const(aBase-BLOCKING8);
     2102          double bValue = bBase[i];
     2103          b0 -= bValue * aBase[0];
     2104          b1 -= bValue * aBase[1];
     2105          b2 -= bValue * aBase[2];
     2106          b3 -= bValue * aBase[3];
     2107          b4 -= bValue * aBase[4];
     2108          b5 -= bValue * aBase[5];
     2109          b6 -= bValue * aBase[6];
     2110          b7 -= bValue * aBase[7];
     2111        }
     2112        bBase2[0] = b0;
     2113        bBase2[1] = b1;
     2114        bBase2[2] = b2;
     2115        bBase2[3] = b3;
     2116        bBase2[4] = b4;
     2117        bBase2[5] = b5;
     2118        bBase2[6] = b6;
     2119        bBase2[7] = b7;
    21442120#endif
    21452121#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 }
    2166 void CoinAbcDtrsm1(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
    2167 {
    2168   assert ((m&(BLOCKING8-1))==0);
     2122        __m256d b0 = _mm256_load_pd(bBase2);
     2123        __m256d b1 = _mm256_load_pd(bBase2 + 4);
     2124        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2125          aBase -= BLOCKING8;
     2126          //coin_prefetch_const(aBase5-BLOCKING8);
     2127          __m256d bb = _mm256_broadcast_sd(bBase + i);
     2128          __m256d a0 = _mm256_load_pd(aBase);
     2129          b0 -= a0 * bb;
     2130          __m256d a1 = _mm256_load_pd(aBase + 4);
     2131          b1 -= a1 * bb;
     2132        }
     2133        //_mm256_store_pd (bBase2, b0);
     2134        *reinterpret_cast< __m256d * >(bBase2) = b0;
     2135        //_mm256_store_pd (bBase2+4, b1);
     2136        *reinterpret_cast< __m256d * >(bBase2 + 4) = b1;
     2137#endif
     2138      }
     2139    }
     2140  }
     2141}
     2142void CoinAbcDtrsm1(int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     2143{
     2144  assert((m & (BLOCKING8 - 1)) == 0);
    21692145  // 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);
     2146  for (int kkk = m - BLOCKING8; kkk >= 0; kkk -= UNROLL_DTRSM * BLOCKING8) {
     2147    int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2148    for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2149      const double *COIN_RESTRICT aBase = a + m * ii + ii * BLOCKING8;
     2150      double *COIN_RESTRICT bBase = b + ii;
     2151      for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2152        bBase[i] *= aBase[i * (BLOCKING8 + 1)];
     2153        for (int k = i - 1; k >= 0; k--) {
     2154          bBase[k] -= bBase[i] * aBase[k + i * BLOCKING8];
     2155        }
     2156      }
     2157      double *COIN_RESTRICT bBase2 = bBase;
     2158      for (int iii = ii; iii > mm; iii -= BLOCKING8) {
     2159        bBase2 -= BLOCKING8;
     2160        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2161          aBase -= BLOCKING8;
     2162          coin_prefetch_const(aBase - BLOCKING8);
     2163          for (int k = BLOCKING8 - 1; k >= 0; k--) {
     2164            bBase2[k] -= bBase[i] * aBase[k];
     2165          }
     2166        }
     2167      }
     2168    }
     2169    dtrsm1(kkk, 0, mm, m, a, b);
    21942170  }
    21952171}
    21962172static 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);
     2173  int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     2174{
     2175  int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2176  assert((last - first) % BLOCKING8 == 0);
     2177  if (last - first > CILK_DTRSM) {
     2178    int mid = ((first + last) >> 4) << 3;
     2179    cilk_spawn dtrsm2(kkk, first, mid, m, a, b);
     2180    dtrsm2(kkk, mid, last, m, a, b);
    22052181    cilk_sync;
    22062182  } 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 }
    2259 void CoinAbcDtrsm2(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
    2260 {
    2261   assert ((m&(BLOCKING8-1))==0);
     2183    for (int iii = last - BLOCKING8; iii >= first; iii -= BLOCKING8) {
     2184      for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2185        const double *COIN_RESTRICT bBase = b + ii;
     2186        double *COIN_RESTRICT bBase2 = b + iii;
     2187        const double *COIN_RESTRICT aBase = a + ii * BLOCKING8 + m * iii + BLOCKING8X8;
     2188        for (int j = BLOCKING8 - 1; j >= 0; j -= 4) {
     2189          double bValue0 = bBase2[j];
     2190          double bValue1 = bBase2[j - 1];
     2191          double bValue2 = bBase2[j - 2];
     2192          double bValue3 = bBase2[j - 3];
     2193          bValue0 -= aBase[-1 * BLOCKING8 + 7] * bBase[7];
     2194          bValue1 -= aBase[-2 * BLOCKING8 + 7] * bBase[7];
     2195          bValue2 -= aBase[-3 * BLOCKING8 + 7] * bBase[7];
     2196          bValue3 -= aBase[-4 * BLOCKING8 + 7] * bBase[7];
     2197          bValue0 -= aBase[-1 * BLOCKING8 + 6] * bBase[6];
     2198          bValue1 -= aBase[-2 * BLOCKING8 + 6] * bBase[6];
     2199          bValue2 -= aBase[-3 * BLOCKING8 + 6] * bBase[6];
     2200          bValue3 -= aBase[-4 * BLOCKING8 + 6] * bBase[6];
     2201          bValue0 -= aBase[-1 * BLOCKING8 + 5] * bBase[5];
     2202          bValue1 -= aBase[-2 * BLOCKING8 + 5] * bBase[5];
     2203          bValue2 -= aBase[-3 * BLOCKING8 + 5] * bBase[5];
     2204          bValue3 -= aBase[-4 * BLOCKING8 + 5] * bBase[5];
     2205          bValue0 -= aBase[-1 * BLOCKING8 + 4] * bBase[4];
     2206          bValue1 -= aBase[-2 * BLOCKING8 + 4] * bBase[4];
     2207          bValue2 -= aBase[-3 * BLOCKING8 + 4] * bBase[4];
     2208          bValue3 -= aBase[-4 * BLOCKING8 + 4] * bBase[4];
     2209          bValue0 -= aBase[-1 * BLOCKING8 + 3] * bBase[3];
     2210          bValue1 -= aBase[-2 * BLOCKING8 + 3] * bBase[3];
     2211          bValue2 -= aBase[-3 * BLOCKING8 + 3] * bBase[3];
     2212          bValue3 -= aBase[-4 * BLOCKING8 + 3] * bBase[3];
     2213          bValue0 -= aBase[-1 * BLOCKING8 + 2] * bBase[2];
     2214          bValue1 -= aBase[-2 * BLOCKING8 + 2] * bBase[2];
     2215          bValue2 -= aBase[-3 * BLOCKING8 + 2] * bBase[2];
     2216          bValue3 -= aBase[-4 * BLOCKING8 + 2] * bBase[2];
     2217          bValue0 -= aBase[-1 * BLOCKING8 + 1] * bBase[1];
     2218          bValue1 -= aBase[-2 * BLOCKING8 + 1] * bBase[1];
     2219          bValue2 -= aBase[-3 * BLOCKING8 + 1] * bBase[1];
     2220          bValue3 -= aBase[-4 * BLOCKING8 + 1] * bBase[1];
     2221          bValue0 -= aBase[-1 * BLOCKING8 + 0] * bBase[0];
     2222          bValue1 -= aBase[-2 * BLOCKING8 + 0] * bBase[0];
     2223          bValue2 -= aBase[-3 * BLOCKING8 + 0] * bBase[0];
     2224          bValue3 -= aBase[-4 * BLOCKING8 + 0] * bBase[0];
     2225          bBase2[j] = bValue0;
     2226          bBase2[j - 1] = bValue1;
     2227          bBase2[j - 2] = bValue2;
     2228          bBase2[j - 3] = bValue3;
     2229          aBase -= 4 * BLOCKING8;
     2230        }
     2231      }
     2232    }
     2233  }
     2234}
     2235void CoinAbcDtrsm2(int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     2236{
     2237  assert((m & (BLOCKING8 - 1)) == 0);
    22622238  // 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);
     2239  for (int kkk = m - BLOCKING8; kkk >= 0; kkk -= UNROLL_DTRSM * BLOCKING8) {
     2240    int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2241    for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2242      double *COIN_RESTRICT bBase = b + ii;
     2243      double *COIN_RESTRICT bBase2 = bBase;
     2244      const double *COIN_RESTRICT aBase = a + m * ii + (ii + BLOCKING8) * BLOCKING8;
     2245      for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2246        aBase -= BLOCKING8;
     2247        for (int k = i + 1; k < BLOCKING8; k++) {
     2248          bBase2[i] -= aBase[k] * bBase2[k];
     2249        }
     2250      }
     2251      for (int iii = ii - BLOCKING8; iii >= mm; iii -= BLOCKING8) {
     2252        bBase2 -= BLOCKING8;
     2253        assert(bBase2 == b + iii);
     2254        aBase -= m * BLOCKING8;
     2255        const double *COIN_RESTRICT aBase2 = aBase + BLOCKING8X8;
     2256        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2257          aBase2 -= BLOCKING8;
     2258          for (int k = BLOCKING8 - 1; k >= 0; k--) {
     2259            bBase2[i] -= aBase2[k] * bBase[k];
     2260          }
     2261        }
     2262      }
     2263    }
     2264    dtrsm2(kkk, 0, mm, m, a, b);
    22892265  }
    22902266}
     
    22922268#define CILK_DTRSM3 32
    22932269static void dtrsm3(int kkk, int first, int last,
    2294                    int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
     2270  int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
    22952271{
    22962272  //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);
     2273  assert((last - first) % BLOCKING8 == 0);
     2274  if (last - first > CILK_DTRSM3) {
     2275    int mid = ((first + last) >> 4) << 3;
     2276    cilk_spawn dtrsm3(kkk, first, mid, m, a, b);
     2277    dtrsm3(kkk, mid, last, m, a, b);
    23022278    cilk_sync;
    23032279  } 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 }
    2356 void CoinAbcDtrsm3(int m, double * COIN_RESTRICT a,double * COIN_RESTRICT b)
    2357 {
    2358   assert ((m&(BLOCKING8-1))==0);
     2280    for (int kk = 0; kk < kkk; kk += BLOCKING8) {
     2281      for (int ii = first; ii < last; ii += BLOCKING8) {
     2282        double *COIN_RESTRICT bBase = b + ii;
     2283        double *COIN_RESTRICT bBase2 = b + kk;
     2284        const double *COIN_RESTRICT aBase = a + ii * m + kk * BLOCKING8;
     2285        for (int j = 0; j < BLOCKING8; j += 4) {
     2286          double bValue0 = bBase[j];
     2287          double bValue1 = bBase[j + 1];
     2288          double bValue2 = bBase[j + 2];
     2289          double bValue3 = bBase[j + 3];
     2290          bValue0 -= aBase[0] * bBase2[0];
     2291          bValue1 -= aBase[1 * BLOCKING8 + 0] * bBase2[0];
     2292          bValue2 -= aBase[2 * BLOCKING8 + 0] * bBase2[0];
     2293          bValue3 -= aBase[3 * BLOCKING8 + 0] * bBase2[0];
     2294          bValue0 -= aBase[1] * bBase2[1];
     2295          bValue1 -= aBase[1 * BLOCKING8 + 1] * bBase2[1];
     2296          bValue2 -= aBase[2 * BLOCKING8 + 1] * bBase2[1];
     2297          bValue3 -= aBase[3 * BLOCKING8 + 1] * bBase2[1];
     2298          bValue0 -= aBase[2] * bBase2[2];
     2299          bValue1 -= aBase[1 * BLOCKING8 + 2] * bBase2[2];
     2300          bValue2 -= aBase[2 * BLOCKING8 + 2] * bBase2[2];
     2301          bValue3 -= aBase[3 * BLOCKING8 + 2] * bBase2[2];
     2302          bValue0 -= aBase[3] * bBase2[3];
     2303          bValue1 -= aBase[1 * BLOCKING8 + 3] * bBase2[3];
     2304          bValue2 -= aBase[2 * BLOCKING8 + 3] * bBase2[3];
     2305          bValue3 -= aBase[3 * BLOCKING8 + 3] * bBase2[3];
     2306          bValue0 -= aBase[4] * bBase2[4];
     2307          bValue1 -= aBase[1 * BLOCKING8 + 4] * bBase2[4];
     2308          bValue2 -= aBase[2 * BLOCKING8 + 4] * bBase2[4];
     2309          bValue3 -= aBase[3 * BLOCKING8 + 4] * bBase2[4];
     2310          bValue0 -= aBase[5] * bBase2[5];
     2311          bValue1 -= aBase[1 * BLOCKING8 + 5] * bBase2[5];
     2312          bValue2 -= aBase[2 * BLOCKING8 + 5] * bBase2[5];
     2313          bValue3 -= aBase[3 * BLOCKING8 + 5] * bBase2[5];
     2314          bValue0 -= aBase[6] * bBase2[6];
     2315          bValue1 -= aBase[1 * BLOCKING8 + 6] * bBase2[6];
     2316          bValue2 -= aBase[2 * BLOCKING8 + 7] * bBase2[7];
     2317          bValue3 -= aBase[3 * BLOCKING8 + 6] * bBase2[6];
     2318          bValue0 -= aBase[7] * bBase2[7];
     2319          bValue1 -= aBase[1 * BLOCKING8 + 7] * bBase2[7];
     2320          bValue2 -= aBase[2 * BLOCKING8 + 6] * bBase2[6];
     2321          bValue3 -= aBase[3 * BLOCKING8 + 7] * bBase2[7];
     2322          bBase[j] = bValue0;
     2323          bBase[j + 1] = bValue1;
     2324          bBase[j + 2] = bValue2;
     2325          bBase[j + 3] = bValue3;
     2326          aBase += 4 * BLOCKING8;
     2327        }
     2328      }
     2329    }
     2330  }
     2331}
     2332void CoinAbcDtrsm3(int m, double *COIN_RESTRICT a, double *COIN_RESTRICT b)
     2333{
     2334  assert((m & (BLOCKING8 - 1)) == 0);
    23592335  // 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 }
    2386 static void
    2387 CoinAbcDlaswp(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 }
    2413 extern 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                          );
    2419 static 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);
     2336  for (int kkk = 0; kkk < m; kkk += UNROLL_DTRSM3 * BLOCKING8) {
     2337    int mm = CoinMin(m, kkk + UNROLL_DTRSM3 * BLOCKING8);
     2338    dtrsm3(kkk, kkk, mm, m, a, b);
     2339    for (int ii = kkk; ii < mm; ii += BLOCKING8) {
     2340      double *COIN_RESTRICT bBase = b + ii;
     2341      for (int kk = kkk; kk < ii; kk += BLOCKING8) {
     2342        double *COIN_RESTRICT bBase2 = b + kk;
     2343        const double *COIN_RESTRICT aBase = a + ii * m + kk * BLOCKING8;
     2344        for (int i = 0; i < BLOCKING8; i++) {
     2345          for (int k = 0; k < BLOCKING8; k++) {
     2346            bBase[i] -= aBase[k] * bBase2[k];
     2347          }
     2348          aBase += BLOCKING8;
     2349        }
     2350      }
     2351      const double *COIN_RESTRICT aBase = a + ii * m + ii * BLOCKING8;
     2352      for (int i = 0; i < BLOCKING8; i++) {
     2353        for (int k = 0; k < i; k++) {
     2354          bBase[i] -= aBase[k] * bBase[k];
     2355        }
     2356        bBase[i] = bBase[i] * aBase[i];
     2357        aBase += BLOCKING8;
     2358      }
     2359    }
     2360  }
     2361}
     2362static void
     2363CoinAbcDlaswp(int n, double *COIN_RESTRICT a, int lda, int start, int end, int *ipiv)
     2364{
     2365  assert((n & (BLOCKING8 - 1)) == 0);
     2366  double *COIN_RESTRICT aThis3 = a;
     2367  for (int j = 0; j < n; j += BLOCKING8) {
     2368    for (int i = start; i < end; i++) {
     2369      int ip = ipiv[i];
     2370      if (ip != i) {
     2371        double *COIN_RESTRICT aTop = aThis3 + j * lda;
     2372        int iBias = ip / BLOCKING8;
     2373        ip -= iBias * BLOCKING8;
     2374        double *COIN_RESTRICT aNotTop = aTop + iBias * BLOCKING8X8;
     2375        iBias = i / BLOCKING8;
     2376        int i2 = i - iBias * BLOCKING8;
     2377        aTop += iBias * BLOCKING8X8;
     2378        for (int k = 0; k < BLOCKING8; k++) {
     2379          double temp = aTop[i2];
     2380          aTop[i2] = aNotTop[ip];
     2381          aNotTop[ip] = temp;
     2382          aTop += BLOCKING8;
     2383          aNotTop += BLOCKING8;
     2384        }
     2385      }
     2386    }
     2387  }
     2388}
     2389extern void CoinAbcDgemm(int m, int n, int k, double *COIN_RESTRICT a, int lda,
     2390  double *COIN_RESTRICT b, double *COIN_RESTRICT c
     2391#if ABC_PARALLEL == 2
     2392  ,
     2393  int parallelMode
     2394#endif
     2395);
     2396static int CoinAbcDgetrf2(int m, int n, double *COIN_RESTRICT a, int *ipiv)
     2397{
     2398  assert((m & (BLOCKING8 - 1)) == 0 && (n & (BLOCKING8 - 1)) == 0);
     2399  assert(n == BLOCKING8);
     2400  int dimension = CoinMin(m, n);
    24242401  /* entry for column j and row i (when multiple of BLOCKING8)
    24252402     is at aBlocked+j*m+i*BLOCKING8
    24262403  */
    2427   assert (dimension==BLOCKING8);
     2404  assert(dimension == BLOCKING8);
    24282405  //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;
     2406  double *COIN_RESTRICT aThis3 = a;
     2407  for (int j = 0; j < dimension; j++) {
     2408    int pivotRow = -1;
     2409    double largest = 0.0;
     2410    double *COIN_RESTRICT aThis2 = aThis3 + j * BLOCKING8;
    24342411    // 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;
     2412    int pivotRow2 = -1;
     2413    for (int i = j; i < BLOCKING8; i++) {
     2414      double value = fabs(aThis2[i]);
     2415      if (value > largest) {
     2416        largest = value;
     2417        pivotRow2 = i;
    24412418      }
    24422419    }
    24432420    // 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;
     2421    int iBias = 0;
     2422    for (int ii = BLOCKING8; ii < m; ii += BLOCKING8) {
     2423      aThis2 += BLOCKING8 * BLOCKING8;
     2424      for (int i = 0; i < BLOCKING8; i++) {
     2425        double value = fabs(aThis2[i]);
     2426        if (value > largest) {
     2427          largest = value;
     2428          pivotRow2 = i;
     2429          iBias = ii;
     2430        }
     2431      }
     2432    }
     2433    pivotRow = pivotRow2 + iBias;
     2434    ipiv[j] = pivotRow;
    24582435    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;
     2436      if (j != pivotRow) {
     2437        double *COIN_RESTRICT aTop = aThis3;
     2438        double *COIN_RESTRICT aNotTop = aThis3 + iBias * BLOCKING8;
     2439        for (int i = 0; i < n; i++) {
     2440          double value = aTop[j];
     2441          aTop[j] = aNotTop[pivotRow2];
     2442          aNotTop[pivotRow2] = value;
     2443          aTop += BLOCKING8;
     2444          aNotTop += BLOCKING8;
     2445        }
     2446      }
     2447      aThis2 = aThis3 + j * BLOCKING8;
    24712448      double pivotMultiplier = 1.0 / aThis2[j];
    24722449      aThis2[j] = pivotMultiplier;
    24732450      // Do L
    24742451      // this block
    2475       for (int i=j+1;i<BLOCKING8;i++)
    2476         aThis2[i] *= pivotMultiplier;
     2452      for (int i = j + 1; i < BLOCKING8; i++)
     2453        aThis2[i] *= pivotMultiplier;
    24772454      // 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;
     2455      for (int ii = BLOCKING8; ii < m; ii += BLOCKING8) {
     2456        aThis2 += BLOCKING8 * BLOCKING8;
     2457        for (int i = 0; i < BLOCKING8; i++)
     2458          aThis2[i] *= pivotMultiplier;
    24822459      }
    24832460      // update rest
    2484       double * COIN_RESTRICT aPut;
    2485       aThis2 = aThis3+j*BLOCKING8;
    2486       aPut = aThis2+BLOCKING8;
    2487       double * COIN_RESTRICT aPut2 = aPut;
     2461      double *COIN_RESTRICT aPut;
     2462      aThis2 = aThis3 + j * BLOCKING8;
     2463      aPut = aThis2 + BLOCKING8;
     2464      double *COIN_RESTRICT aPut2 = aPut;
    24882465      // 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;
     2466      for (int i = j + 1; i < BLOCKING8; i++) {
     2467        double value = aPut[j];
     2468        for (int k = j + 1; k < BLOCKING8; k++)
     2469          aPut[k] -= value * aThis2[k];
     2470        aPut += BLOCKING8;
    24942471      }
    24952472      // 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         }
     2473      for (int ii = BLOCKING8; ii < m; ii += BLOCKING8) {
     2474        aThis2 += BLOCKING8 * BLOCKING8;
     2475        aPut = aThis2 + BLOCKING8;
     2476        double *COIN_RESTRICT aPut2a = aPut2;
     2477        for (int i = j + 1; i < BLOCKING8; i++) {
     2478          double value = aPut2a[j];
     2479          for (int k = 0; k < BLOCKING8; k++)
     2480            aPut[k] -= value * aThis2[k];
     2481          aPut += BLOCKING8;
     2482          aPut2a += BLOCKING8;
     2483        }
    25072484      }
    25082485    } else {
     
    25132490}
    25142491
    2515 int
    2516 CoinAbcDgetrf(int m, int n, double * COIN_RESTRICT a, int lda, int * ipiv
    2517 #if ABC_PARALLEL==2
    2518                           ,int parallelMode
     2492int CoinAbcDgetrf(int m, int n, double *COIN_RESTRICT a, int lda, int *ipiv
     2493#if ABC_PARALLEL == 2
     2494  ,
     2495  int parallelMode
    25192496#endif
    25202497)
    25212498{
    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);
     2499  assert(m == n);
     2500  assert((m & (BLOCKING8 - 1)) == 0 && (n & (BLOCKING8 - 1)) == 0);
     2501  if (m < BLOCKING8) {
     2502    return CoinAbcDgetrf2(m, n, a, ipiv);
    25262503  } 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);
     2504    for (int j = 0; j < n; j += BLOCKING8) {
     2505      int start = j;
     2506      int newSize = CoinMin(BLOCKING8, n - j);
     2507      int end = j + newSize;
     2508      int returnCode = CoinAbcDgetrf2(m - start, newSize, a + (start * lda + start * BLOCKING8),
     2509        ipiv + start);
    25332510      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         }
     2511        // adjust
     2512        for (int k = start; k < end; k++)
     2513          ipiv[k] += start;
     2514        // swap 0<start
     2515        CoinAbcDlaswp(start, a, lda, start, end, ipiv);
     2516        if (end < n) {
     2517          // swap >=end
     2518          CoinAbcDlaswp(n - end, a + end * lda, lda, start, end, ipiv);
     2519          CoinAbcDtrsmFactor(newSize, n - end, a + (start * lda + start * BLOCKING8), lda);
     2520          CoinAbcDgemm(n - end, n - end, newSize,
     2521            a + start * lda + end * BLOCKING8, lda,
     2522            a + end * lda + start * BLOCKING8, a + end * lda + end * BLOCKING8
     2523#if ABC_PARALLEL == 2
     2524            ,
     2525            parallelMode
     2526#endif
     2527          );
     2528        }
    25512529      } else {
    2552         return returnCode;
     2530        return returnCode;
    25532531      }
    25542532    }
     
    25562534  return 0;
    25572535}
    2558 void
    2559 CoinAbcDgetrs(char trans,int m, double * COIN_RESTRICT a, double * COIN_RESTRICT work)
    2560 {
    2561   assert ((m&(BLOCKING8-1))==0);
    2562   if (trans=='N') {
     2536void CoinAbcDgetrs(char trans, int m, double *COIN_RESTRICT a, double *COIN_RESTRICT work)
     2537{
     2538  assert((m & (BLOCKING8 - 1)) == 0);
     2539  if (trans == 'N') {
    25632540    //CoinAbcDlaswp1(work,m,ipiv);
    2564     CoinAbcDtrsm0(m,a,work);
    2565     CoinAbcDtrsm1(m,a,work);
     2541    CoinAbcDtrsm0(m, a, work);
     2542    CoinAbcDtrsm1(m, a, work);
    25662543  } else {
    2567     assert (trans=='T');
    2568     CoinAbcDtrsm3(m,a,work);
    2569     CoinAbcDtrsm2(m,a,work);
     2544    assert(trans == 'T');
     2545    CoinAbcDtrsm3(m, a, work);
     2546    CoinAbcDtrsm2(m, a, work);
    25702547    //CoinAbcDlaswp1Backwards(work,m,ipiv);
    25712548  }
     
    25792556   3 Left Upper Transpose NonUnit
    25802557*/
    2581 static 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);
     2558static void CoinAbcDtrsmFactor(int m, int n, long double *COIN_RESTRICT a, int lda)
     2559{
     2560  assert((m & (BLOCKING8 - 1)) == 0 && (n & (BLOCKING8 - 1)) == 0);
     2561  assert(m == BLOCKING8);
    25852562  // 0 Left Lower NoTranspose Unit
    25862563  /* entry for column j and row i (when multiple of BLOCKING8)
    25872564     is at aBlocked+j*m+i*BLOCKING8
    25882565  */
    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++) {
     2566  long double *COIN_RESTRICT aBase2 = a;
     2567  long double *COIN_RESTRICT bBase2 = aBase2 + lda * BLOCKING8;
     2568  for (int jj = 0; jj < n; jj += BLOCKING8) {
     2569    long double *COIN_RESTRICT bBase = bBase2;
     2570    for (int j = jj; j < jj + BLOCKING8; j++) {
    25942571#if 0
    25952572      long double * COIN_RESTRICT aBase = aBase2;
     
    26052582#else
    26062583      // unroll - stay in registers - don't test for zeros
    2607       assert (BLOCKING8==8);
     2584      assert(BLOCKING8 == 8);
    26082585      long double bValue0 = bBase[0];
    26092586      long double bValue1 = bBase[1];
     
    26142591      long double bValue6 = bBase[6];
    26152592      long double bValue7 = bBase[7];
    2616       bValue1-=bValue0*a[1+0*BLOCKING8];
     2593      bValue1 -= bValue0 * a[1 + 0 * BLOCKING8];
    26172594      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];
     2595      bValue2 -= bValue0 * a[2 + 0 * BLOCKING8];
     2596      bValue3 -= bValue0 * a[3 + 0 * BLOCKING8];
     2597      bValue4 -= bValue0 * a[4 + 0 * BLOCKING8];
     2598      bValue5 -= bValue0 * a[5 + 0 * BLOCKING8];
     2599      bValue6 -= bValue0 * a[6 + 0 * BLOCKING8];
     2600      bValue7 -= bValue0 * a[7 + 0 * BLOCKING8];
     2601      bValue2 -= bValue1 * a[2 + 1 * BLOCKING8];
    26252602      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];
     2603      bValue3 -= bValue1 * a[3 + 1 * BLOCKING8];
     2604      bValue4 -= bValue1 * a[4 + 1 * BLOCKING8];
     2605      bValue5 -= bValue1 * a[5 + 1 * BLOCKING8];
     2606      bValue6 -= bValue1 * a[6 + 1 * BLOCKING8];
     2607      bValue7 -= bValue1 * a[7 + 1 * BLOCKING8];
     2608      bValue3 -= bValue2 * a[3 + 2 * BLOCKING8];
    26322609      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];
     2610      bValue4 -= bValue2 * a[4 + 2 * BLOCKING8];
     2611      bValue5 -= bValue2 * a[5 + 2 * BLOCKING8];
     2612      bValue6 -= bValue2 * a[6 + 2 * BLOCKING8];
     2613      bValue7 -= bValue2 * a[7 + 2 * BLOCKING8];
     2614      bValue4 -= bValue3 * a[4 + 3 * BLOCKING8];
    26382615      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];
     2616      bValue5 -= bValue3 * a[5 + 3 * BLOCKING8];
     2617      bValue6 -= bValue3 * a[6 + 3 * BLOCKING8];
     2618      bValue7 -= bValue3 * a[7 + 3 * BLOCKING8];
     2619      bValue5 -= bValue4 * a[5 + 4 * BLOCKING8];
    26432620      bBase[5] = bValue5;
    2644       bValue6-=bValue4*a[6+4*BLOCKING8];
    2645       bValue7-=bValue4*a[7+4*BLOCKING8];
    2646       bValue6-=bValue5*a[6+5*BLOCKING8];
     2621      bValue6 -= bValue4 * a[6 + 4 * BLOCKING8];
     2622      bValue7 -= bValue4 * a[7 + 4 * BLOCKING8];
     2623      bValue6 -= bValue5 * a[6 + 5 * BLOCKING8];
    26472624      bBase[6] = bValue6;
    2648       bValue7-=bValue5*a[7+5*BLOCKING8];
    2649       bValue7-=bValue6*a[7+6*BLOCKING8];
     2625      bValue7 -= bValue5 * a[7 + 5 * BLOCKING8];
     2626      bValue7 -= bValue6 * a[7 + 6 * BLOCKING8];
    26502627      bBase[7] = bValue7;
    26512628#endif
    26522629      bBase += BLOCKING8;
    26532630    }
    2654     bBase2 +=lda*BLOCKING8;
     2631    bBase2 += lda * BLOCKING8;
    26552632  }
    26562633}
     
    26582635#define CILK_DTRSM 32
    26592636static 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);
     2637  int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     2638{
     2639  int mm = CoinMin(kkk + UNROLL_DTRSM * BLOCKING8, m);
     2640  assert((last - first) % BLOCKING8 == 0);
     2641  if (last - first > CILK_DTRSM) {
     2642    int mid = ((first + last) >> 4) << 3;
     2643    cilk_spawn dtrsm0(kkk, first, mid, m, a, b);
     2644    dtrsm0(kkk, mid, last, m, a, b);
    26682645    cilk_sync;
    26692646  } 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) {
     2647    const long double *COIN_RESTRICT aBaseA = a + UNROLL_DTRSM * BLOCKING8X8 + kkk * BLOCKING8;
     2648    aBaseA += (first - mm) * BLOCKING8 - BLOCKING8X8;
     2649    aBaseA += m * kkk;
     2650    for (int ii = first; ii < last; ii += BLOCKING8) {
    26742651      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
     2652      const long double *COIN_RESTRICT aBaseB = aBaseA;
     2653      long double *COIN_RESTRICT bBaseI = b + ii;
     2654      for (int kk = kkk; kk < mm; kk += BLOCKING8) {
     2655        long double *COIN_RESTRICT bBase = b + kk;
     2656        const long double *COIN_RESTRICT aBase2 = aBaseB; //a+UNROLL_DTRSM*BLOCKING8X8+m*kk+kkk*BLOCKING8;
     2657        //aBase2 += (ii-mm)*BLOCKING8;
     2658        //assert (aBase2==aBaseB);
     2659        aBaseB += m * BLOCKING8;
     2660#if AVX2 != 2
    26842661#define ALTERNATE_INNER
    26852662#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         }
     2663        for (int k = 0; k < BLOCKING8; k++) {
     2664          //coin_prefetch_const(aBase2+BLOCKING8);
     2665          for (int i = 0; i < BLOCKING8; i++) {
     2666            bBaseI[i] -= bBase[k] * aBase2[i];
     2667          }
     2668          aBase2 += BLOCKING8;
     2669        }
    26932670#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;
     2671        long double b0 = bBaseI[0];
     2672        long double b1 = bBaseI[1];
     2673        long double b2 = bBaseI[2];
     2674        long double b3 = bBaseI[3];
     2675        long double b4 = bBaseI[4];
     2676        long double b5 = bBaseI[5];
     2677        long double b6 = bBaseI[6];
     2678        long double b7 = bBaseI[7];
     2679        for (int k = 0; k < BLOCKING8; k++) {
     2680          //coin_prefetch_const(aBase2+BLOCKING8);
     2681          long double bValue = bBase[k];
     2682          b0 -= bValue * aBase2[0];
     2683          b1 -= bValue * aBase2[1];
     2684          b2 -= bValue * aBase2[2];
     2685          b3 -= bValue * aBase2[3];
     2686          b4 -= bValue * aBase2[4];
     2687          b5 -= bValue * aBase2[5];
     2688          b6 -= bValue * aBase2[6];
     2689          b7 -= bValue * aBase2[7];
     2690          aBase2 += BLOCKING8;
     2691        }
     2692        bBaseI[0] = b0;
     2693        bBaseI[1] = b1;
     2694        bBaseI[2] = b2;
     2695        bBaseI[3] = b3;
     2696        bBaseI[4] = b4;
     2697        bBaseI[5] = b5;
     2698        bBaseI[6] = b6;
     2699        bBaseI[7] = b7;
    27232700#endif
    27242701#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 }
    2745 void CoinAbcDtrsm0(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
    2746 {
    2747   assert ((m&(BLOCKING8-1))==0);
     2702        __m256d b0 = _mm256_load_pd(bBaseI);
     2703        __m256d b1 = _mm256_load_pd(bBaseI + 4);
     2704        for (int k = 0; k < BLOCKING8; k++) {
     2705          //coin_prefetch_const(aBase2+BLOCKING8);
     2706          __m256d bb = _mm256_broadcast_sd(bBase + k);
     2707          __m256d a0 = _mm256_load_pd(aBase2);
     2708          b0 -= a0 * bb;
     2709          __m256d a1 = _mm256_load_pd(aBase2 + 4);
     2710          b1 -= a1 * bb;
     2711          aBase2 += BLOCKING8;
     2712        }
     2713        //_mm256_store_pd ((bBaseI), (b0));
     2714        *reinterpret_cast< __m256d * >(bBaseI) = b0;
     2715        //_mm256_store_pd (bBaseI+4, b1);
     2716        *reinterpret_cast< __m256d * >(bBaseI + 4) = b1;
     2717#endif
     2718      }
     2719    }
     2720  }
     2721}
     2722void CoinAbcDtrsm0(int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     2723{
     2724  assert((m & (BLOCKING8 - 1)) == 0);
    27482725  // 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);
     2726  for (int kkk = 0; kkk < m; kkk += UNROLL_DTRSM * BLOCKING8) {
     2727    int mm = CoinMin(m, kkk + UNROLL_DTRSM * BLOCKING8);
     2728    for (int kk = kkk; kk < mm; kk += BLOCKING8) {
     2729      const long double *COIN_RESTRICT aBase2 = a + kk * (m + BLOCKING8);
     2730      long double *COIN_RESTRICT bBase = b + kk;
     2731      for (int k = 0; k < BLOCKING8; k++) {
     2732        for (int i = k + 1; i < BLOCKING8; i++) {
     2733          bBase[i] -= bBase[k] * aBase2[i];
     2734        }
     2735        aBase2 += BLOCKING8;
     2736      }
     2737      for (int ii = kk + BLOCKING8; ii < mm; ii += BLOCKING8) {
     2738        long double *COIN_RESTRICT bBaseI = b + ii;
     2739        for (int k = 0; k < BLOCKING8; k++) {
     2740          //coin_prefetch_const(aBase2+BLOCKING8);
     2741          for (int i = 0; i < BLOCKING8; i++) {
     2742            bBaseI[i] -= bBase[k] * aBase2[i];
     2743          }
     2744          aBase2 += BLOCKING8;
     2745        }
     2746      }
     2747    }
     2748    dtrsm0(kkk, mm, m, m, a, b);
    27722749  }
    27732750}
    27742751static 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);
     2752  int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     2753{
     2754  int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2755  assert((last - first) % BLOCKING8 == 0);
     2756  if (last - first > CILK_DTRSM) {
     2757    int mid = ((first + last) >> 4) << 3;
     2758    cilk_spawn dtrsm1(kkk, first, mid, m, a, b);
     2759    dtrsm1(kkk, mid, last, m, a, b);
    27832760    cilk_sync;
    27842761  } 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
     2762    for (int iii = last - BLOCKING8; iii >= first; iii -= BLOCKING8) {
     2763      long double *COIN_RESTRICT bBase2 = b + iii;
     2764      const long double *COIN_RESTRICT aBaseA = a + BLOCKING8X8 + BLOCKING8 * iii;
     2765      for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2766        long double *COIN_RESTRICT bBase = b + ii;
     2767        const long double *COIN_RESTRICT aBase = aBaseA + m * ii;
     2768#if AVX2 != 2
    27932769#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         }
     2770        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2771          aBase -= BLOCKING8;
     2772          //coin_prefetch_const(aBase-BLOCKING8);
     2773          for (int k = BLOCKING8 - 1; k >= 0; k--) {
     2774            bBase2[k] -= bBase[i] * aBase[k];
     2775          }
     2776        }
    28012777#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;
     2778        long double b0 = bBase2[0];
     2779        long double b1 = bBase2[1];
     2780        long double b2 = bBase2[2];
     2781        long double b3 = bBase2[3];
     2782        long double b4 = bBase2[4];
     2783        long double b5 = bBase2[5];
     2784        long double b6 = bBase2[6];
     2785        long double b7 = bBase2[7];
     2786        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2787          aBase -= BLOCKING8;
     2788          //coin_prefetch_const(aBase-BLOCKING8);
     2789          long double bValue = bBase[i];
     2790          b0 -= bValue * aBase[0];
     2791          b1 -= bValue * aBase[1];
     2792          b2 -= bValue * aBase[2];
     2793          b3 -= bValue * aBase[3];
     2794          b4 -= bValue * aBase[4];
     2795          b5 -= bValue * aBase[5];
     2796          b6 -= bValue * aBase[6];
     2797          b7 -= bValue * aBase[7];
     2798        }
     2799        bBase2[0] = b0;
     2800        bBase2[1] = b1;
     2801        bBase2[2] = b2;
     2802        bBase2[3] = b3;
     2803        bBase2[4] = b4;
     2804        bBase2[5] = b5;
     2805        bBase2[6] = b6;
     2806        bBase2[7] = b7;
    28312807#endif
    28322808#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 }
    2853 void CoinAbcDtrsm1(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
    2854 {
    2855   assert ((m&(BLOCKING8-1))==0);
     2809        __m256d b0 = _mm256_load_pd(bBase2);
     2810        __m256d b1 = _mm256_load_pd(bBase2 + 4);
     2811        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2812          aBase -= BLOCKING8;
     2813          //coin_prefetch_const(aBase5-BLOCKING8);
     2814          __m256d bb = _mm256_broadcast_sd(bBase + i);
     2815          __m256d a0 = _mm256_load_pd(aBase);
     2816          b0 -= a0 * bb;
     2817          __m256d a1 = _mm256_load_pd(aBase + 4);
     2818          b1 -= a1 * bb;
     2819        }
     2820        //_mm256_store_pd (bBase2, b0);
     2821        *reinterpret_cast< __m256d * >(bBase2) = b0;
     2822        //_mm256_store_pd (bBase2+4, b1);
     2823        *reinterpret_cast< __m256d * >(bBase2 + 4) = b1;
     2824#endif
     2825      }
     2826    }
     2827  }
     2828}
     2829void CoinAbcDtrsm1(int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     2830{
     2831  assert((m & (BLOCKING8 - 1)) == 0);
    28562832  // 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);
     2833  for (int kkk = m - BLOCKING8; kkk >= 0; kkk -= UNROLL_DTRSM * BLOCKING8) {
     2834    int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2835    for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2836      const long double *COIN_RESTRICT aBase = a + m * ii + ii * BLOCKING8;
     2837      long double *COIN_RESTRICT bBase = b + ii;
     2838      for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2839        bBase[i] *= aBase[i * (BLOCKING8 + 1)];
     2840        for (int k = i - 1; k >= 0; k--) {
     2841          bBase[k] -= bBase[i] * aBase[k + i * BLOCKING8];
     2842        }
     2843      }
     2844      long double *COIN_RESTRICT bBase2 = bBase;
     2845      for (int iii = ii; iii > mm; iii -= BLOCKING8) {
     2846        bBase2 -= BLOCKING8;
     2847        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2848          aBase -= BLOCKING8;
     2849          coin_prefetch_const(aBase - BLOCKING8);
     2850          for (int k = BLOCKING8 - 1; k >= 0; k--) {
     2851            bBase2[k] -= bBase[i] * aBase[k];
     2852          }
     2853        }
     2854      }
     2855    }
     2856    dtrsm1(kkk, 0, mm, m, a, b);
    28812857  }
    28822858}
    28832859static 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);
     2860  int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     2861{
     2862  int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2863  assert((last - first) % BLOCKING8 == 0);
     2864  if (last - first > CILK_DTRSM) {
     2865    int mid = ((first + last) >> 4) << 3;
     2866    cilk_spawn dtrsm2(kkk, first, mid, m, a, b);
     2867    dtrsm2(kkk, mid, last, m, a, b);
    28922868    cilk_sync;
    28932869  } 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 }
    2946 void CoinAbcDtrsm2(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
    2947 {
    2948   assert ((m&(BLOCKING8-1))==0);
     2870    for (int iii = last - BLOCKING8; iii >= first; iii -= BLOCKING8) {
     2871      for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2872        const long double *COIN_RESTRICT bBase = b + ii;
     2873        long double *COIN_RESTRICT bBase2 = b + iii;
     2874        const long double *COIN_RESTRICT aBase = a + ii * BLOCKING8 + m * iii + BLOCKING8X8;
     2875        for (int j = BLOCKING8 - 1; j >= 0; j -= 4) {
     2876          long double bValue0 = bBase2[j];
     2877          long double bValue1 = bBase2[j - 1];
     2878          long double bValue2 = bBase2[j - 2];
     2879          long double bValue3 = bBase2[j - 3];
     2880          bValue0 -= aBase[-1 * BLOCKING8 + 7] * bBase[7];
     2881          bValue1 -= aBase[-2 * BLOCKING8 + 7] * bBase[7];
     2882          bValue2 -= aBase[-3 * BLOCKING8 + 7] * bBase[7];
     2883          bValue3 -= aBase[-4 * BLOCKING8 + 7] * bBase[7];
     2884          bValue0 -= aBase[-1 * BLOCKING8 + 6] * bBase[6];
     2885          bValue1 -= aBase[-2 * BLOCKING8 + 6] * bBase[6];
     2886          bValue2 -= aBase[-3 * BLOCKING8 + 6] * bBase[6];
     2887          bValue3 -= aBase[-4 * BLOCKING8 + 6] * bBase[6];
     2888          bValue0 -= aBase[-1 * BLOCKING8 + 5] * bBase[5];
     2889          bValue1 -= aBase[-2 * BLOCKING8 + 5] * bBase[5];
     2890          bValue2 -= aBase[-3 * BLOCKING8 + 5] * bBase[5];
     2891          bValue3 -= aBase[-4 * BLOCKING8 + 5] * bBase[5];
     2892          bValue0 -= aBase[-1 * BLOCKING8 + 4] * bBase[4];
     2893          bValue1 -= aBase[-2 * BLOCKING8 + 4] * bBase[4];
     2894          bValue2 -= aBase[-3 * BLOCKING8 + 4] * bBase[4];
     2895          bValue3 -= aBase[-4 * BLOCKING8 + 4] * bBase[4];
     2896          bValue0 -= aBase[-1 * BLOCKING8 + 3] * bBase[3];
     2897          bValue1 -= aBase[-2 * BLOCKING8 + 3] * bBase[3];
     2898          bValue2 -= aBase[-3 * BLOCKING8 + 3] * bBase[3];
     2899          bValue3 -= aBase[-4 * BLOCKING8 + 3] * bBase[3];
     2900          bValue0 -= aBase[-1 * BLOCKING8 + 2] * bBase[2];
     2901          bValue1 -= aBase[-2 * BLOCKING8 + 2] * bBase[2];
     2902          bValue2 -= aBase[-3 * BLOCKING8 + 2] * bBase[2];
     2903          bValue3 -= aBase[-4 * BLOCKING8 + 2] * bBase[2];
     2904          bValue0 -= aBase[-1 * BLOCKING8 + 1] * bBase[1];
     2905          bValue1 -= aBase[-2 * BLOCKING8 + 1] * bBase[1];
     2906          bValue2 -= aBase[-3 * BLOCKING8 + 1] * bBase[1];
     2907          bValue3 -= aBase[-4 * BLOCKING8 + 1] * bBase[1];
     2908          bValue0 -= aBase[-1 * BLOCKING8 + 0] * bBase[0];
     2909          bValue1 -= aBase[-2 * BLOCKING8 + 0] * bBase[0];
     2910          bValue2 -= aBase[-3 * BLOCKING8 + 0] * bBase[0];
     2911          bValue3 -= aBase[-4 * BLOCKING8 + 0] * bBase[0];
     2912          bBase2[j] = bValue0;
     2913          bBase2[j - 1] = bValue1;
     2914          bBase2[j - 2] = bValue2;
     2915          bBase2[j - 3] = bValue3;
     2916          aBase -= 4 * BLOCKING8;
     2917        }
     2918      }
     2919    }
     2920  }
     2921}
     2922void CoinAbcDtrsm2(int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     2923{
     2924  assert((m & (BLOCKING8 - 1)) == 0);
    29492925  // 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);
     2926  for (int kkk = m - BLOCKING8; kkk >= 0; kkk -= UNROLL_DTRSM * BLOCKING8) {
     2927    int mm = CoinMax(0, kkk - (UNROLL_DTRSM - 1) * BLOCKING8);
     2928    for (int ii = kkk; ii >= mm; ii -= BLOCKING8) {
     2929      long double *COIN_RESTRICT bBase = b + ii;
     2930      long double *COIN_RESTRICT bBase2 = bBase;
     2931      const long double *COIN_RESTRICT aBase = a + m * ii + (ii + BLOCKING8) * BLOCKING8;
     2932      for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2933        aBase -= BLOCKING8;
     2934        for (int k = i + 1; k < BLOCKING8; k++) {
     2935          bBase2[i] -= aBase[k] * bBase2[k];
     2936        }
     2937      }
     2938      for (int iii = ii - BLOCKING8; iii >= mm; iii -= BLOCKING8) {
     2939        bBase2 -= BLOCKING8;
     2940        assert(bBase2 == b + iii);
     2941        aBase -= m * BLOCKING8;
     2942        const long double *COIN_RESTRICT aBase2 = aBase + BLOCKING8X8;
     2943        for (int i = BLOCKING8 - 1; i >= 0; i--) {
     2944          aBase2 -= BLOCKING8;
     2945          for (int k = BLOCKING8 - 1; k >= 0; k--) {
     2946            bBase2[i] -= aBase2[k] * bBase[k];
     2947          }
     2948        }
     2949      }
     2950    }
     2951    dtrsm2(kkk, 0, mm, m, a, b);
    29762952  }
    29772953}
     
    29792955#define CILK_DTRSM3 32
    29802956static void dtrsm3(int kkk, int first, int last,
    2981                    int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
     2957  int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
    29822958{
    29832959  //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);
     2960  assert((last - first) % BLOCKING8 == 0);
     2961  if (last - first > CILK_DTRSM3) {
     2962    int mid = ((first + last) >> 4) << 3;
     2963    cilk_spawn dtrsm3(kkk, first, mid, m, a, b);
     2964    dtrsm3(kkk, mid, last, m, a, b);
    29892965    cilk_sync;
    29902966  } 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 }
    3043 void CoinAbcDtrsm3(int m, long double * COIN_RESTRICT a,long double * COIN_RESTRICT b)
    3044 {
    3045   assert ((m&(BLOCKING8-1))==0);
     2967    for (int kk = 0; kk < kkk; kk += BLOCKING8) {
     2968      for (int ii = first; ii < last; ii += BLOCKING8) {
     2969        long double *COIN_RESTRICT bBase = b + ii;
     2970        long double *COIN_RESTRICT bBase2 = b + kk;
     2971        const long double *COIN_RESTRICT aBase = a + ii * m + kk * BLOCKING8;
     2972        for (int j = 0; j < BLOCKING8; j += 4) {
     2973          long double bValue0 = bBase[j];
     2974          long double bValue1 = bBase[j + 1];
     2975          long double bValue2 = bBase[j + 2];
     2976          long double bValue3 = bBase[j + 3];
     2977          bValue0 -= aBase[0] * bBase2[0];
     2978          bValue1 -= aBase[1 * BLOCKING8 + 0] * bBase2[0];
     2979          bValue2 -= aBase[2 * BLOCKING8 + 0] * bBase2[0];
     2980          bValue3 -= aBase[3 * BLOCKING8 + 0] * bBase2[0];
     2981          bValue0 -= aBase[1] * bBase2[1];
     2982          bValue1 -= aBase[1 * BLOCKING8 + 1] * bBase2[1];
     2983          bValue2 -= aBase[2 * BLOCKING8 + 1] * bBase2[1];
     2984          bValue3 -= aBase[3 * BLOCKING8 + 1] * bBase2[1];
     2985          bValue0 -= aBase[2] * bBase2[2];
     2986          bValue1 -= aBase[1 * BLOCKING8 + 2] * bBase2[2];
     2987          bValue2 -= aBase[2 * BLOCKING8 + 2] * bBase2[2];
     2988          bValue3 -= aBase[3 * BLOCKING8 + 2] * bBase2[2];
     2989          bValue0 -= aBase[3] * bBase2[3];
     2990          bValue1 -= aBase[1 * BLOCKING8 + 3] * bBase2[3];
     2991          bValue2 -= aBase[2 * BLOCKING8 + 3] * bBase2[3];
     2992          bValue3 -= aBase[3 * BLOCKING8 + 3] * bBase2[3];
     2993          bValue0 -= aBase[4] * bBase2[4];
     2994          bValue1 -= aBase[1 * BLOCKING8 + 4] * bBase2[4];
     2995          bValue2 -= aBase[2 * BLOCKING8 + 4] * bBase2[4];
     2996          bValue3 -= aBase[3 * BLOCKING8 + 4] * bBase2[4];
     2997          bValue0 -= aBase[5] * bBase2[5];
     2998          bValue1 -= aBase[1 * BLOCKING8 + 5] * bBase2[5];
     2999          bValue2 -= aBase[2 * BLOCKING8 + 5] * bBase2[5];
     3000          bValue3 -= aBase[3 * BLOCKING8 + 5] * bBase2[5];
     3001          bValue0 -= aBase[6] * bBase2[6];
     3002          bValue1 -= aBase[1 * BLOCKING8 + 6] * bBase2[6];
     3003          bValue2 -= aBase[2 * BLOCKING8 + 7] * bBase2[7];
     3004          bValue3 -= aBase[3 * BLOCKING8 + 6] * bBase2[6];
     3005          bValue0 -= aBase[7] * bBase2[7];
     3006          bValue1 -= aBase[1 * BLOCKING8 + 7] * bBase2[7];
     3007          bValue2 -= aBase[2 * BLOCKING8 + 6] * bBase2[6];
     3008          bValue3 -= aBase[3 * BLOCKING8 + 7] * bBase2[7];
     3009          bBase[j] = bValue0;
     3010          bBase[j + 1] = bValue1;
     3011          bBase[j + 2] = bValue2;
     3012          bBase[j + 3] = bValue3;
     3013          aBase += 4 * BLOCKING8;
     3014        }
     3015      }
     3016    }
     3017  }
     3018}
     3019void CoinAbcDtrsm3(int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT b)
     3020{
     3021  assert((m & (BLOCKING8 - 1)) == 0);
    30463022  // 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 }
    3073 static void
    3074 CoinAbcDlaswp(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 }
    3100 extern 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                          );
    3106 static 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);
     3023  for (int kkk = 0; kkk < m; kkk += UNROLL_DTRSM3 * BLOCKING8) {
     3024    int mm = CoinMin(m, kkk + UNROLL_DTRSM3 * BLOCKING8);
     3025    dtrsm3(kkk, kkk, mm, m, a, b);
     3026    for (int ii = kkk; ii < mm; ii += BLOCKING8) {
     3027      long double *COIN_RESTRICT bBase = b + ii;
     3028      for (int kk = kkk; kk < ii; kk += BLOCKING8) {
     3029        long double *COIN_RESTRICT bBase2 = b + kk;
     3030        const long double *COIN_RESTRICT aBase = a + ii * m + kk * BLOCKING8;
     3031        for (int i = 0; i < BLOCKING8; i++) {
     3032          for (int k = 0; k < BLOCKING8; k++) {
     3033            bBase[i] -= aBase[k] * bBase2[k];
     3034          }
     3035          aBase += BLOCKING8;
     3036        }
     3037      }
     3038      const long double *COIN_RESTRICT aBase = a + ii * m + ii * BLOCKING8;
     3039      for (int i = 0; i < BLOCKING8; i++) {
     3040        for (int k = 0; k < i; k++) {
     3041          bBase[i] -= aBase[k] * bBase[k];
     3042        }
     3043        bBase[i] = bBase[i] * aBase[i];
     3044        aBase += BLOCKING8;
     3045      }
     3046    }
     3047  }
     3048}
     3049static void
     3050CoinAbcDlaswp(int n, long double *COIN_RESTRICT a, int lda, int start, int end, int *ipiv)
     3051{
     3052  assert((n & (BLOCKING8 - 1)) == 0);
     3053  long double *COIN_RESTRICT aThis3 = a;
     3054  for (int j = 0; j < n; j += BLOCKING8) {
     3055    for (int i = start; i < end; i++) {
     3056      int ip = ipiv[i];
     3057      if (ip != i) {
     3058        long double *COIN_RESTRICT aTop = aThis3 + j * lda;
     3059        int iBias = ip / BLOCKING8;
     3060        ip -= iBias * BLOCKING8;
     3061        long double *COIN_RESTRICT aNotTop = aTop + iBias * BLOCKING8X8;
     3062        iBias = i / BLOCKING8;
     3063        int i2 = i - iBias * BLOCKING8;
     3064        aTop += iBias * BLOCKING8X8;
     3065        for (int k = 0; k < BLOCKING8; k++) {
     3066          long double temp = aTop[i2];
     3067          aTop[i2] = aNotTop[ip];
     3068          aNotTop[ip] = temp;
     3069          aTop += BLOCKING8;
     3070          aNotTop += BLOCKING8;
     3071        }
     3072      }
     3073    }
     3074  }
     3075}
     3076extern void CoinAbcDgemm(int m, int n, int k, long double *COIN_RESTRICT a, int lda,
     3077  long double *COIN_RESTRICT b, long double *COIN_RESTRICT c
     3078#if ABC_PARALLEL == 2
     3079  ,
     3080  int parallelMode
     3081#endif
     3082);
     3083static int CoinAbcDgetrf2(int m, int n, long double *COIN_RESTRICT a, int *ipiv)
     3084{
     3085  assert((m & (BLOCKING8 - 1)) == 0 && (n & (BLOCKING8 - 1)) == 0);
     3086  assert(n == BLOCKING8);
     3087  int dimension = CoinMin(m, n);
    31113088  /* entry for column j and row i (when multiple of BLOCKING8)
    31123089     is at aBlocked+j*m+i*BLOCKING8
    31133090  */
    3114   assert (dimension==BLOCKING8);
     3091  assert(dimension == BLOCKING8);
    31153092  //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;
     3093  long double *COIN_RESTRICT aThis3 = a;
     3094  for (int j = 0; j < dimension; j++) {
     3095    int pivotRow = -1;
     3096    long double largest = 0.0;
     3097    long double *COIN_RESTRICT aThis2 = aThis3 + j * BLOCKING8;
    31213098    // 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;
     3099    int pivotRow2 = -1;
     3100    for (int i = j; i < BLOCKING8; i++) {
     3101      long double value = fabsl(aThis2[i]);
     3102      if (value > largest) {
     3103        largest = value;
     3104        pivotRow2 = i;
    31283105      }
    31293106    }
    31303107    // 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;
     3108    int iBias = 0;
     3109    for (int ii = BLOCKING8; ii < m; ii += BLOCKING8) {
     3110      aThis2 += BLOCKING8 * BLOCKING8;
     3111      for (int i = 0; i < BLOCKING8; i++) {
     3112        long double value = fabsl(aThis2[i]);
     3113        if (value > largest) {
     3114          largest = value;
     3115          pivotRow2 = i;
     3116          iBias = ii;
     3117        }
     3118      }
     3119    }
     3120    pivotRow = pivotRow2 + iBias;
     3121    ipiv[j] = pivotRow;
    31453122    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;
     3123      if (j != pivotRow) {
     3124        long double *COIN_RESTRICT aTop = aThis3;
     3125        long double *COIN_RESTRICT aNotTop = aThis3 + iBias * BLOCKING8;
     3126        for (int i = 0; i < n; i++) {
     3127          long double value = aTop[j];
     3128          aTop[j] = aNotTop[pivotRow2];
     3129          aNotTop[pivotRow2] = value;
     3130          aTop += BLOCKING8;
     3131          aNotTop += BLOCKING8;
     3132        }
     3133      }
     3134      aThis2 = aThis3 + j * BLOCKING8;
    31583135      long double pivotMultiplier = 1.0 / aThis2[j];
    31593136      aThis2[j] = pivotMultiplier;
    31603137      // Do L
    31613138      // this block
    3162       for (int i=j+1;i<BLOCKING8;i++)
    3163         aThis2[i] *= pivotMultiplier;
     3139      for (int i = j + 1; i < BLOCKING8; i++)
     3140        aThis2[i] *= pivotMultiplier;
    31643141      // 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;
     3142      for (int ii = BLOCKING8; ii < m; ii += BLOCKING8) {
     3143        aThis2 += BLOCKING8 * BLOCKING8;
     3144        for (int i = 0; i < BLOCKING8; i++)
     3145          aThis2[i] *= pivotMultiplier;
    31693146      }
    31703147      // update rest
    3171       long double * COIN_RESTRICT aPut;
    3172       aThis2 = aThis3+j*BLOCKING8;
    3173       aPut = aThis2+BLOCKING8;
    3174       long double * COIN_RESTRICT aPut2 = aPut;
     3148      long double *COIN_RESTRICT aPut;
     3149      aThis2 = aThis3 + j * BLOCKING8;
     3150      aPut = aThis2 + BLOCKING8;
     3151      long double *COIN_RESTRICT aPut2 = aPut;
    31753152      // 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;
     3153      for (int i = j + 1; i < BLOCKING8; i++) {
     3154        long double value = aPut[j];
     3155        for (int k = j + 1; k < BLOCKING8; k++)
     3156          aPut[k] -= value * aThis2[k];
     3157        aPut += BLOCKING8;
    31813158      }
    31823159      // 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         }
     3160      for (int ii = BLOCKING8; ii < m; ii += BLOCKING8) {
     3161        aThis2 += BLOCKING8 * BLOCKING8;
     3162        aPut = aThis2 + BLOCKING8;
     3163        long double *COIN_RESTRICT aPut2a = aPut2;
     3164        for (int i = j + 1; i < BLOCKING8; i++) {
     3165          long double value = aPut2a[j];
     3166          for (int k = 0; k < BLOCKING8; k++)
     3167            aPut[k] -= value * aThis2[k];
     3168          aPut += BLOCKING8;
     3169          aPut2a += BLOCKING8;
     3170        }
    31943171      }
    31953172    } else {
     
    32003177}
    32013178
    3202 int
    3203 CoinAbcDgetrf(int m, int n, long double * COIN_RESTRICT a, int lda, int * ipiv
    3204 #if ABC_PARALLEL==2
    3205                           ,int parallelMode
     3179int CoinAbcDgetrf(int m, int n, long double *COIN_RESTRICT a, int lda, int *ipiv
     3180#if ABC_PARALLEL == 2
     3181  ,
     3182  int parallelMode
    32063183#endif
    32073184)
    32083185{
    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);
     3186  assert(m == n);
     3187  assert((m & (BLOCKING8 - 1)) == 0 && (n & (BLOCKING8 - 1)) == 0);
     3188  if (m < BLOCKING8) {
     3189    return CoinAbcDgetrf2(m, n, a, ipiv);
    32133190  } 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);
     3191    for (int j = 0; j < n; j += BLOCKING8) {
     3192      int start = j;
     3193      int newSize = CoinMin(BLOCKING8, n - j);
     3194      int end = j + newSize;
     3195      int returnCode = CoinAbcDgetrf2(m - start, newSize, a + (start * lda + start * BLOCKING8),
     3196        ipiv + start);
    32203197      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         }
     3198        // adjust
     3199        for (int k = start; k < end; k++)
     3200          ipiv[k] += start;
     3201        // swap 0<start
     3202        CoinAbcDlaswp(start, a, lda, start, end, ipiv);
     3203        if (end < n) {
     3204          // swap >=end
     3205          CoinAbcDlaswp(n - end, a + end * lda, lda, start, end, ipiv);
     3206          CoinAbcDtrsmFactor(newSize, n - end, a + (start * lda + start * BLOCKING8), lda);
     3207          CoinAbcDgemm(n - end, n - end, newSize,
     3208            a + start * lda + end * BLOCKING8, lda,
     3209            a + end * lda + start * BLOCKING8, a + end * lda + end * BLOCKING8
     3210#if ABC_PARALLEL == 2
     3211            ,
     3212            parallelMode
     3213#endif
     3214          );
     3215        }
    32383216      } else {
    3239         return returnCode;
     3217        return returnCode;
    32403218      }
    32413219    }
     
    32433221  return 0;
    32443222}
    3245 void
    3246 CoinAbcDgetrs(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') {
     3223void CoinAbcDgetrs(char trans, int m, long double *COIN_RESTRICT a, long double *COIN_RESTRICT work)
     3224{
     3225  assert((m & (BLOCKING8 - 1)) == 0);
     3226  if (trans == 'N') {
    32503227    //CoinAbcDlaswp1(work,m,ipiv);
    3251     CoinAbcDtrsm0(m,a,work);
    3252     CoinAbcDtrsm1(m,a,work);
     3228    CoinAbcDtrsm0(m, a, work);
     3229    CoinAbcDtrsm1(m, a, work);
    32533230  } else {
    3254     assert (trans=='T');
    3255     CoinAbcDtrsm3(m,a,work);
    3256     CoinAbcDtrsm2(m,a,work);
     3231    assert(trans == 'T');
     3232    CoinAbcDtrsm3(m, a, work);
     3233    CoinAbcDtrsm2(m, a, work);
    32573234    //CoinAbcDlaswp1Backwards(work,m,ipiv);
    32583235  }
    32593236}
    32603237#endif
     3238
     3239/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     3240*/
Note: See TracChangeset for help on using the changeset viewer.