Changeset 2259 for trunk


Ignore:
Timestamp:
Mar 27, 2017 4:37:20 AM (3 years ago)
Author:
forrest
Message:

for avx stuff

Location:
trunk/Clp/src
Files:
13 edited

Legend:

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

    r2235 r2259  
    710710               model->setMoreSpecialOptions(value);
    711711               break;
     712          case CLP_PARAM_INT_VECTOR_MODE:
     713               model->setVectorMode(value);
     714               break;
    712715#ifndef COIN_HAS_CBC
    713716#ifdef CBC_THREAD
     
    751754     case CLP_PARAM_INT_MORESPECIALOPTIONS:
    752755          value = model->moreSpecialOptions();
     756          break;
     757     case CLP_PARAM_INT_VECTOR_MODE:
     758          value = model->vectorMode();
    753759          break;
    754760#ifndef COIN_HAS_CBC
     
    39253931It is possible you can get same effect by using example driver4.cpp."
    39263932     );
     3933#endif
     3934#ifdef COIN_AVX2
     3935     parameters[numberParameters++] =
     3936          CbcOrClpParam("vector!Mode", "Try and use vector instructions",
     3937                        0, 1, CLP_PARAM_INT_VECTOR_MODE);
     3938     parameters[numberParameters-1].setLonghelp
     3939     (
     3940          "At present only for Intel architectures - but could be extended.  \
     3941Uses avx2 or avx512 instructions. Uses different storage for matrix - can be \
     3942of benefit without instruction set on some problems."
     3943     );
     3944     parameters[numberParameters-1].setIntValue(0);
     3945#endif
     3946#ifdef COIN_HAS_CBC
    39273947     parameters[numberParameters++] =
    39283948          CbcOrClpParam("Vnd!VariableNeighborhoodSearch", "Whether to try Variable Neighborhood Search",
  • trunk/Clp/src/CbcOrClpParam.hpp

    r2176 r2259  
    105105     CLP_PARAM_INT_MORESPECIALOPTIONS,
    106106     CLP_PARAM_INT_DECOMPOSE_BLOCKS,
     107     CLP_PARAM_INT_VECTOR_MODE,
    107108
    108109     CBC_PARAM_INT_STRONGBRANCHING = 151,
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r2257 r2259  
    751751          return;
    752752     }
     753     if (columnCopy_)
     754       factor *= 0.7;
    753755     if (numberInRowArray > factor * numberRows || !rowCopy) {
    754756          // do by column
     
    10321034                              double tentativeTheta = 1.0e15;
    10331035                              double upperTheta = 1.0e31;
    1034                               double bestPossible = 0.0;
    10351036                              int addSequence = model->numberColumns();
    10361037                              const unsigned char * COIN_RESTRICT statusArray = model->statusArray() + addSequence;
     
    10491050                                             value = oldValue - tentativeTheta * alpha;
    10501051                                             if (value < dualT) {
    1051                                                   bestPossible = CoinMax(bestPossible, alpha);
    10521052                                                  value = oldValue - upperTheta * alpha;
    10531053                                                  if (value < dualT && alpha >= acceptablePivot) {
     
    10711071                                             model->djRegion(1),
    10721072                                             upperTheta,
    1073                                              bestPossible,
    10741073                                             acceptablePivot,
    10751074                                             model->currentDualTolerance(),
     
    10771076                                             zeroTolerance);
    10781077                              model->spareDoubleArray_[0] = upperTheta;
    1079                               model->spareDoubleArray_[1] = bestPossible;
    10801078                              spareArray->setNumElements(numberRemaining);
    10811079#if 0
     
    10851083                              spareArray->checkClean();
    10861084                              columnArray->checkClean();
    1087                               printf("col %d spare %d upper %g best %g\n",
     1085                              printf("col %d spare %d upper %g\n",
    10881086                                     columnArray->getNumElements(),
    1089                                      spareArray->getNumElements(),upperTheta,bestPossible);
     1087                                     spareArray->getNumElements(),upperTheta);
    10901088#endif
    10911089                              // signal partially done
     
    11091107                    //xA++;
    11101108               } else {
     1109                 if ((model->moreSpecialOptions()&8)!=0 &&
     1110                     model->algorithm()<0) {
     1111                   columnCopy_->transposeTimes(model, pi, columnArray,
     1112                     model->rowArray(3),rowArray);
     1113                   model->spareIntArray_[0]=-2;
     1114                 } else {
    11111115                    columnCopy_->transposeTimes(model, pi, columnArray);
     1116                 }
    11121117                    numberNonZero = columnArray->getNumElements();
    11131118                    //xB++;
     
    11501155                    //xC++;
    11511156               } else {
     1157                 if ((model->moreSpecialOptions()&8)!=0 &&
     1158                     model->algorithm()<0) {
     1159                   columnCopy_->transposeTimes(model, pi, columnArray,
     1160                     model->rowArray(3),rowArray);
     1161                   model->spareIntArray_[0]=-2;
     1162                 } else {
    11521163                    columnCopy_->transposeTimes(model, pi, columnArray);
     1164                 }
    11531165                    numberNonZero = columnArray->getNumElements();
    11541166                    //xD++;
     
    17381750  int numberRemaining = 0;
    17391751  double upperTheta = info.upperTheta;
    1740   double bestPossible = info.bestPossible;
    17411752  int first = info.startColumn;
    17421753  int last = first+info.numberToDo;
     
    17761787          double value = oldValue - tentativeTheta * alpha;
    17771788          if (value < dualT) {
    1778             bestPossible = CoinMax(bestPossible, alpha);
    17791789            value = oldValue - upperTheta * alpha;
    17801790            if (value < dualT && alpha >= acceptablePivot) {
     
    17921802  info.numberAdded=numberNonZero;
    17931803  info.numberRemaining=numberRemaining;
    1794   info.bestPossible=bestPossible;
    17951804  info.upperTheta=upperTheta;
    17961805}
     
    18071816          const double * COIN_RESTRICT reducedCost,
    18081817          double & upperThetaP,
    1809           double & bestPossibleP,
    18101818          double acceptablePivot,
    18111819          double dualTolerance,
     
    18151823     int numberRemaining = numberRemainingP;
    18161824     double upperTheta = upperThetaP;
    1817      double bestPossible = bestPossibleP;
    18181825     int numberNonZero = 0;
    18191826     // get matrix data pointers
     
    18291836     for (int i=0;i<numberThreads;i++) {
    18301837       info[i].upperTheta = upperThetaP;
    1831        info[i].bestPossible = bestPossibleP;
    18321838       info[i].acceptablePivot = acceptablePivot;
    18331839       info[i].reducedCost = const_cast<double *>(reducedCost);
     
    18541860       numberNonZero += info[i].numberAdded;
    18551861       numberRemaining += info[i].numberRemaining;
    1856        bestPossible = CoinMax(bestPossible,static_cast<double>(info[i].bestPossible));
    18571862       upperTheta = CoinMin(upperTheta,static_cast<double>(info[i].upperTheta));
    18581863     }
     
    19741979                         double value = oldValue - tentativeTheta * alpha;
    19751980                         if (value < dualT) {
    1976                               bestPossible = CoinMax(bestPossible, alpha);
    19771981                              value = oldValue - upperTheta * alpha;
    19781982                              if (value < dualT && alpha >= acceptablePivot) {
     
    19931997     numberRemainingP = numberRemaining;
    19941998     upperThetaP = upperTheta;
    1995      bestPossibleP = bestPossible;
    19961999     return numberNonZero;
    19972000}
     
    58685871{
    58695872     delete columnCopy_;
    5870      if ((model->moreSpecialOptions()&16777216)!=0) {
     5873     if (model->vectorMode()) {
    58715874       flags_ |= 16;
    58725875       // go to exact devex (unless full steepest)
     
    58765879         pricing->setMode(0);
    58775880     }
    5878      if ((flags_ & 16) != 0) {
     5881     if ((flags_ & 16) != 0 && model->numberRows()>200 &&
     5882         model->numberColumns()>500) {
    58795883          columnCopy_ = new ClpPackedMatrix3(model, matrix_);
    58805884          flags_ |= 8;
     
    59045908               columnCopy_->sortBlocks(model);
    59055909          }
     5910#ifndef NDEBUG
     5911          columnCopy_->checkBlocks(model);
     5912#endif
    59065913     }
    59075914}
     
    62056212                       int * spareIndex, const double * arrayTemp,
    62066213                       const int * indexTemp, int numberIn,
    6207                        int offset, double acceptablePivot, double * bestPossiblePtr,
     6214                       int offset, double acceptablePivot,
    62086215                       double * upperThetaPtr, int * posFreePtr, double * freePivotPtr)
    62096216{
     
    62116218     int i;
    62126219     int numberRemaining = 0;
    6213      double bestPossible = 0.0;
    62146220     double upperTheta = 1.0e31;
    62156221     double freePivot = acceptablePivot;
     
    62336239          case ClpSimplex::isFree:
    62346240          case ClpSimplex::superBasic:
    6235                bestPossible = CoinMax(bestPossible, fabs(alpha));
    62366241               oldValue = reducedCost[iSequence];
    62376242               // If free has to be very large - should come in via dualRow
     
    62616266               //assert (oldValue<=dualTolerance*1.0001);
    62626267               if (value > dualTolerance) {
    6263                     bestPossible = CoinMax(bestPossible, -alpha);
    62646268                    value = oldValue - upperTheta * alpha;
    62656269                    if (value > dualTolerance && -alpha >= acceptablePivot)
     
    62756279               //assert (oldValue>=-dualTolerance*1.0001);
    62766280               if (value < -dualTolerance) {
    6277                     bestPossible = CoinMax(bestPossible, alpha);
    62786281                    value = oldValue - upperTheta * alpha;
    62796282                    if (value < -dualTolerance && alpha >= acceptablePivot)
     
    62866289          }
    62876290     }
    6288      *bestPossiblePtr = bestPossible;
    62896291     *upperThetaPtr = upperTheta;
    62906292     *freePivotPtr = freePivot;
     
    64626464                                          info->spareIndex, (const double *)info->arrayTemp,
    64636465                                          (const int *) info->indexTemp, *(info->numberInPtr),
    6464                                           info->offset, info->acceptablePivot, info->bestPossiblePtr,
     6466                                          info->offset, info->acceptablePivot,
    64656467                                          info->upperThetaPtr, info->posFreePtr, info->freePivotPtr);
    64666468     return NULL;
     
    64806482     bool dualColumn = model->spareIntArray_[0] == 1;
    64816483     double acceptablePivot = model->spareDoubleArray_[0];
    6482      double bestPossible = 0.0;
    64836484     double upperTheta = 1.0e31;
    64846485     double freePivot = acceptablePivot;
     
    65366537               case ClpSimplex::isFree:
    65376538               case ClpSimplex::superBasic:
    6538                     bestPossible = CoinMax(bestPossible, fabs(alpha));
    65396539                    oldValue = reducedCost[iRow];
    65406540                    // If free has to be very large - should come in via dualRow
     
    65646564                    //assert (oldValue<=dualTolerance*1.0001);
    65656565                    if (value > dualTolerance) {
    6566                          bestPossible = CoinMax(bestPossible, -alpha);
    65676566                         value = oldValue - upperTheta * alpha;
    65686567                         if (value > dualTolerance && -alpha >= acceptablePivot)
     
    65786577                    //assert (oldValue>=-dualTolerance*1.0001);
    65796578                    if (value < -dualTolerance) {
    6580                          bestPossible = CoinMax(bestPossible, alpha);
    65816579                         value = oldValue - upperTheta * alpha;
    65826580                         if (value < -dualTolerance && alpha >= acceptablePivot)
     
    66556653                                      arrayTemp, indexTemp,
    66566654                                      iwork[0], offset3, acceptablePivot,
    6657                                       &dwork[0], &dwork[1], &iwork[2],
     6655                                      &dwork[1], &iwork[2],
    66586656                                      &dwork[2]);
    66596657               int number = iwork[0];
     
    66766674               }
    66776675               upperTheta =  CoinMin(dwork[1], upperTheta);
    6678                bestPossible = CoinMax(dwork[0], bestPossible);
    66796676               for (i = 0; i < number; i++) {
    66806677                    // double value = arrayTemp[i];
     
    66986695               infoPtr->offset = offset;
    66996696               infoPtr->acceptablePivot = acceptablePivot;
    6700                infoPtr->bestPossiblePtr = &dwork[0];
    67016697               infoPtr->upperThetaPtr = &dwork[1];
    67026698               infoPtr->posFreePtr = &iwork[2];
     
    67456741               }
    67466742               upperTheta =  CoinMin(dwork[1], upperTheta);
    6747                bestPossible = CoinMax(dwork[0], bestPossible);
    67486743          }
    67496744          double * arrayTemp = array + offset;
     
    67616756     if (dualColumn) {
    67626757          model->spareDoubleArray_[0] = upperTheta;
    6763           model->spareDoubleArray_[1] = bestPossible;
    67646758          // and theta and alpha and sequence
    67656759          if (posFree < 0) {
     
    68216815  element_(NULL),
    68226816  temporary_(NULL),
    6823   block_(NULL)
    6824 {
    6825 }
     6817  block_(NULL),
     6818  ifActive_(0)
     6819{
     6820}
     6821#ifdef _MSC_VER
     6822#include <intrin.h>
     6823#else
     6824#include <immintrin.h>
     6825//#include <fmaintrin.h>
     6826#endif
    68266827/* Constructor from copy. */
    68276828ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model, const CoinPackedMatrix * columnCopy)
     
    68396840  element_(NULL),
    68406841  temporary_(NULL),
    6841   block_(NULL)
    6842 {
     6842  block_(NULL),
     6843  ifActive_(0)
     6844{
     6845  //#undef COIN_AVX2
     6846  //#define COIN_AVX2 8
     6847  //#define NO_AVX_HARDWARE
    68436848#ifndef COIN_AVX2
    68446849#define COIN_AVX2 4
    68456850#else
     6851#if COIN_AVX2==4
     6852#ifndef NO_AVX_HARDWARE
    68466853#define HASWELL
     6854#endif
     6855#elif COIN_AVX2==8
     6856#ifndef NO_AVX_HARDWARE
     6857#define SKYLAKE
     6858#endif
     6859#else
     6860  error
     6861#endif
    68476862#endif
    68486863#if COIN_AVX2 == 1
     
    68536868#define COIN_AVX2_SHIFT 2
    68546869#elif COIN_AVX2 == 8
    6855   error at present;
     6870#define COIN_AVX2_SHIFT 3
    68566871#else
    68576872  error;
     
    71237138  element_(NULL),
    71247139  temporary_(NULL),
    7125   block_(NULL)
     7140  block_(NULL),
     7141  ifActive_(rhs.ifActive_)
    71267142{
    71277143  if (rhs.numberBlocks_) {
     
    71657181    numberElements_ = rhs.numberElements_;
    71667182    maxBlockSize_ = rhs.maxBlockSize_;
     7183    ifActive_ = rhs.ifActive_;
    71677184    if (rhs.numberBlocks_) {
    71687185      block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
     
    71947211ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
    71957212{
     7213  ifActive_=1;
    71967214  int * lookup = column_ + numberColumnsWithGaps_;
    71977215  for (int iBlock = 0; iBlock < numberBlocks_+1; iBlock++) {
     
    72627280    lastPrice = 0;
    72637281    while (lastPrice <= firstNotPrice) {
    7264       // find first lower
     7282      // find first upper
     7283      int iColumn = numberInBlock;
     7284      for (; lastPrice <= firstNotPrice; lastPrice++) {
     7285        iColumn = column[lastPrice];
     7286        if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound)
     7287          break;
     7288      }
     7289      // find last upper
     7290      int jColumn = -1;
     7291      for (; firstNotPrice > lastPrice; firstNotPrice--) {
     7292        jColumn = column[firstNotPrice];
     7293        if (model->getColumnStatus(jColumn) != ClpSimplex::atUpperBound)
     7294          break;
     7295      }
     7296      if (firstNotPrice > lastPrice) {
     7297        assert (column[lastPrice] == iColumn);
     7298        assert (column[firstNotPrice] == jColumn);
     7299        // need to swap
     7300        column[firstNotPrice] = iColumn;
     7301        lookup[iColumn] = firstNotPrice;
     7302        column[lastPrice] = jColumn;
     7303        lookup[jColumn] = lastPrice;
     7304        int startBit = roundDown(lastPrice);
     7305        CoinBigIndex offset =startBit*nel + (lastPrice-startBit);
     7306        double * elementA = element + offset;
     7307        int * rowA = row + offset;
     7308        startBit = roundDown(firstNotPrice);
     7309        offset =startBit*nel + (firstNotPrice-startBit);
     7310        double * elementB = element + offset;
     7311        int * rowB = row + offset;
     7312        for (int i = 0; i < nel*COIN_AVX2; i+=COIN_AVX2) {
     7313          int temp = rowA[i];
     7314          double tempE = elementA[i];
     7315          rowA[i] = rowB[i];
     7316          elementA[i] = elementB[i];
     7317          rowB[i] = temp;
     7318          elementB[i] = tempE;
     7319        }
     7320        firstNotPrice--;
     7321        lastPrice++;
     7322      } else if (lastPrice == firstNotPrice) {
     7323        // make sure correct side
     7324        iColumn = column[lastPrice];
     7325        if (model->getColumnStatus(iColumn) != ClpSimplex::atUpperBound)
     7326          lastPrice++;
     7327        break;
     7328      }
     7329    }
     7330    block->firstAtUpper_ = lastPrice;
     7331    // now lower
     7332    firstNotPrice = lastPrice - 1;
     7333    lastPrice = 0;
     7334    while (lastPrice <= firstNotPrice) {
     7335      // find first basic or fixed
    72657336      int iColumn = numberInBlock;
    72667337      for (; lastPrice <= firstNotPrice; lastPrice++) {
     
    72697340          break;
    72707341      }
    7271       // find last lower
     7342      // find last non basic or fixed
    72727343      int jColumn = -1;
    72737344      for (; firstNotPrice > lastPrice; firstNotPrice--) {
     
    73117382    }
    73127383    block->firstAtLower_ = lastPrice;
    7313     // now upper
    7314     firstNotPrice = lastPrice - 1;
    7315     lastPrice = 0;
    7316     while (lastPrice <= firstNotPrice) {
    7317       // find first basic or fixed
    7318       int iColumn = numberInBlock;
    7319       for (; lastPrice <= firstNotPrice; lastPrice++) {
    7320         iColumn = column[lastPrice];
    7321         if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound)
    7322           break;
    7323       }
    7324       // find last non basic or fixed
    7325       int jColumn = -1;
    7326       for (; firstNotPrice > lastPrice; firstNotPrice--) {
    7327         jColumn = column[firstNotPrice];
    7328         if (model->getColumnStatus(jColumn) != ClpSimplex::atUpperBound)
    7329           break;
    7330       }
    7331       if (firstNotPrice > lastPrice) {
    7332         assert (column[lastPrice] == iColumn);
    7333         assert (column[firstNotPrice] == jColumn);
    7334         // need to swap
    7335         column[firstNotPrice] = iColumn;
    7336         lookup[iColumn] = firstNotPrice;
    7337         column[lastPrice] = jColumn;
    7338         lookup[jColumn] = lastPrice;
    7339         int startBit = roundDown(lastPrice);
    7340         CoinBigIndex offset =startBit*nel + (lastPrice-startBit);
    7341         double * elementA = element + offset;
    7342         int * rowA = row + offset;
    7343         startBit = roundDown(firstNotPrice);
    7344         offset =startBit*nel + (firstNotPrice-startBit);
    7345         double * elementB = element + offset;
    7346         int * rowB = row + offset;
    7347         for (int i = 0; i < nel*COIN_AVX2; i+=COIN_AVX2) {
    7348           int temp = rowA[i];
    7349           double tempE = elementA[i];
    7350           rowA[i] = rowB[i];
    7351           elementA[i] = elementB[i];
    7352           rowB[i] = temp;
    7353           elementB[i] = tempE;
    7354         }
    7355         firstNotPrice--;
    7356         lastPrice++;
    7357       } else if (lastPrice == firstNotPrice) {
    7358         // make sure correct side
    7359         iColumn = column[lastPrice];
    7360         if (model->getColumnStatus(iColumn) != ClpSimplex::atUpperBound)
    7361           lastPrice++;
    7362         break;
    7363       }
    7364     }
    7365     block->firstAtUpper_ = lastPrice;
    73667384#ifndef NDEBUG
    73677385    // check
     
    73727390              model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
    73737391      assert (lookup[iColumn] == i);
    7374       if (i<block->firstAtUpper_) {
     7392      if (i<block->firstAtLower_) {
    73757393        assert (model->getColumnStatus(iColumn) == ClpSimplex::isFree ||
    73767394                model->getColumnStatus(iColumn) == ClpSimplex::superBasic);
    7377       } else if (i<block->firstAtLower_) {
     7395      } else if (i<block->firstAtUpper_) {
     7396        assert (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound);
     7397      } else {
    73787398        assert (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound);
    7379       } else {
    7380         assert (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound);
    73817399      }
    73827400    }
     
    74297447                            int iColumn)
    74307448{
     7449  if (!ifActive_)
     7450    return;
    74317451  int * lookup = column_ + numberColumnsWithGaps_;
    74327452  // position in block
     
    74637483  int from,to;
    74647484  if (kA<block->firstBasic_) {
    7465     if (kA>=block->firstAtLower_) {
     7485    if (kA>=block->firstAtUpper_) {
    74667486      from=2;
    7467     } else if (kA>=block->firstAtUpper_) {
     7487    } else if (kA>=block->firstAtLower_) {
    74687488      from=1;
    74697489    } else {
     
    74767496      model->getColumnStatus(iColumn) == ClpSimplex::isFixed) {
    74777497    to=3;
     7498  } else if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
     7499    to=2;
    74787500  } else if (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound) {
    7479     to=2;
    7480   } else if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
    74817501    to=1;
    74827502  } else {
    74837503    to=0;
    74847504  }
    7485   int * statusCounts = (&block->firstAtUpper_)-1;
     7505  int * statusCounts = (&block->firstAtLower_)-1;
    74867506  if (from<to) {
    74877507    while (from<to) {
     
    75117531    assert (lookup[iColumn] == i);
    75127532    if (model->algorithm()>0) {
    7513       if (i<block->firstAtUpper_) {
     7533      if (i<block->firstAtLower_) {
    75147534        assert (model->getColumnStatus(iColumn) == ClpSimplex::isFree ||
    75157535                model->getColumnStatus(iColumn) == ClpSimplex::superBasic);
    7516       } else if (i<block->firstAtLower_) {
     7536      } else if (i<block->firstAtUpper_) {
     7537        assert (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound);
     7538      } else {
    75177539        assert (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound);
    7518       } else {
    7519         assert (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound);
    75207540      }
    75217541    }
     
    75317551#endif
    75327552}
     7553#ifndef NDEBUG
     7554/* Debug - check blocks */
     7555void
     7556ClpPackedMatrix3::checkBlocks(const ClpSimplex * model)
     7557{
     7558  if (!ifActive_)
     7559    return;
     7560  for (int iBlock=0;iBlock<numberBlocks_+1;iBlock++) {
     7561    blockStruct * block = block_ + iBlock;
     7562    int * column = column_ + block->startIndices_;
     7563    for (int i=0;i<block->firstAtLower_;i++) {
     7564      int iSequence=column[i];
     7565      assert (model->getColumnStatus(iSequence) == ClpSimplex::isFree ||
     7566              model->getColumnStatus(iSequence) != ClpSimplex::superBasic);
     7567    }
     7568    for (int i=block->firstAtLower_;i<block->firstAtUpper_;i++) {
     7569      int iSequence=column[i];
     7570      assert (model->getColumnStatus(iSequence) == ClpSimplex::atLowerBound);
     7571    }
     7572    for (int i=block->firstAtUpper_;i<block->firstBasic_;i++) {
     7573      int iSequence=column[i];
     7574      assert (model->getColumnStatus(iSequence) == ClpSimplex::atUpperBound);
     7575    }
     7576    for (int i=block->firstBasic_;i<block->numberInBlock_;i++) {
     7577      int iSequence=column[i];
     7578      assert (model->getColumnStatus(iSequence) == ClpSimplex::basic ||
     7579              model->getColumnStatus(iSequence) == ClpSimplex::isFixed);
     7580    }
     7581  }
     7582}
     7583#endif
    75337584/* Return <code>x * -1 * A in <code>z</code>.
    75347585   Note - x packed and z will be packed mode
     
    75917642    int nBlock=numberPrice>>COIN_AVX2_SHIFT;
    75927643    numberPrice -= roundDown(numberPrice);
     7644#if defined(HASWELL) || defined(SKYLAKE)
     7645    double * newValues = roundUpDouble((array+numberNonZero));
     7646#ifdef HASWELL
     7647    assert(COIN_AVX2==4);
     7648    for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7649      __m256d arrayX = _mm256_setzero_pd();
     7650      for (int j=0;j<nel;j++) {
     7651        __m128i rows = _mm_load_si128((const __m128i *)row); // was loadu
     7652        __m256d elements = _mm256_load_pd(element);
     7653        __m256d pis = _mm256_i32gather_pd(pi,rows,8);
     7654        arrayX = _mm256_fmadd_pd(pis,elements,arrayX);
     7655        row+=4;
     7656        element+=4;
     7657      }
     7658      _mm256_store_pd(newValues,arrayX);
     7659#else
     7660    assert(COIN_AVX2==8);
     7661    for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7662      __m512d arrayX = _mm512_setzero_pd();
     7663      for (int j=0;j<nel;j++) {
     7664        __m256i rows = _mm256_load_si256((const __m256i *)row);
     7665        __m512d elements = _mm512_load_pd(element);
     7666        __m512d pis = _mm512_i32gather_pd(rows,pi,8);
     7667        arrayX = _mm512_fmadd_pd(pis,elements,arrayX);
     7668        row+=4;
     7669        element+=4;
     7670      }
     7671      _mm512_store_pd(newValues,arrayX);
     7672#endif
     7673      newValues += COIN_AVX2;
     7674      //row += (nel-1)*COIN_AVX2;
     7675      //element += (nel-1)*COIN_AVX2;
     7676      assert (row == row_ + block->startElements_ + nel*COIN_AVX2*(jBlock+1));
     7677    }
     7678    int n = nBlock*COIN_AVX2;
     7679    int nSave=newValues-array;
     7680    newValues = roundUpDouble((array+numberNonZero));
     7681    for (int j=0;j<n;j++) {
     7682      double value=newValues[j];
     7683      if (fabs(value) > zeroTolerance) {
     7684        array[numberNonZero] = value;
     7685        index[numberNonZero++] = *column;
     7686      }
     7687      column++;
     7688    }
     7689    for (int j=numberNonZero;j<nSave;j++)
     7690      array[j]=0.0;
     7691#else
    75937692    for (int jBlock=0;jBlock<nBlock;jBlock++) {
    75947693      for (int j=0; j < COIN_AVX2; j++) {
     
    76177716#endif
    76187717    }
     7718#endif
    76197719    // last lot
    76207720    for (int j=0; j < numberPrice; j++) {
     
    76397739  }
    76407740  output->setNumElements(numberNonZero);
     7741}
     7742/* Return <code>x * -1 * A in <code>z</code>.
     7743   Note - x packed and z will be packed mode
     7744   Squashes small elements and knows about ClpSimplex
     7745   - does dualColumn0 */
     7746void
     7747ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
     7748                                 const double * COIN_RESTRICT pi,
     7749                                 CoinIndexedVector * output,
     7750                                 CoinIndexedVector * candidate,
     7751                                 const CoinIndexedVector * rowArray) const
     7752{
     7753  int numberNonZero = 0;
     7754  int * index = output->getIndices();
     7755  double * array = output->denseVector();
     7756  double zeroTolerance = model->zeroTolerance();
     7757  int numberColumns=model->numberColumns();
     7758  const unsigned char * COIN_RESTRICT statusArray=model->statusArray() +
     7759    numberColumns;
     7760  int numberInRowArray = rowArray->getNumElements();
     7761  const int * COIN_RESTRICT whichRow = rowArray->getIndices();
     7762  const double * COIN_RESTRICT piOld = rowArray->denseVector();
     7763  int * COIN_RESTRICT spareIndex = candidate->getIndices();
     7764  double * COIN_RESTRICT spareArray = candidate->denseVector();
     7765  const double * COIN_RESTRICT reducedCost=model->djRegion(0);
     7766  double multiplier[] = { -1.0, 1.0};
     7767  double dualT = - model->currentDualTolerance();
     7768  double acceptablePivot = model->spareDoubleArray_[0];
     7769  double tentativeTheta = 1.0e15;
     7770  double upperTheta = 1.0e31;
     7771  int numberRemaining=0;
     7772  // dualColumn0 for slacks
     7773  for (int i = 0; i < numberInRowArray; i++) {
     7774    int iSequence = whichRow[i];
     7775    int iStatus = (statusArray[iSequence] & 3) - 1;
     7776    if (iStatus) {
     7777      double mult = multiplier[iStatus-1];
     7778      double alpha = piOld[i] * mult;
     7779      double oldValue;
     7780      double value;
     7781      if (alpha > 0.0) {
     7782        oldValue = reducedCost[iSequence] * mult;
     7783        value = oldValue - tentativeTheta * alpha;
     7784        if (value < dualT) {
     7785          value = oldValue - upperTheta * alpha;
     7786          if (value < dualT && alpha >= acceptablePivot) {
     7787            upperTheta = (oldValue - dualT) / alpha;
     7788          }
     7789          // add to list
     7790          spareArray[numberRemaining] = alpha * mult;
     7791          spareIndex[numberRemaining++] = iSequence + numberColumns;
     7792        }
     7793      }
     7794    }
     7795  }
     7796  statusArray -= numberColumns;
     7797  reducedCost -= numberColumns;
     7798  double value = 0.0;
     7799  CoinBigIndex j;
     7800  int numberOdd = block_->startIndices_;
     7801  if (numberOdd) {
     7802    // A) as probably long may be worth unrolling
     7803    CoinBigIndex end = start_[1];
     7804    for (j = start_[0]; j < end; j++) {
     7805      int iRow = row_[j];
     7806      value += pi[iRow] * element_[j];
     7807    }
     7808    int iColumn;
     7809    // int jColumn=column_[0];
     7810   
     7811    for (iColumn = 0; iColumn < numberOdd - 1; iColumn++) {
     7812      CoinBigIndex start = end;
     7813      end = start_[iColumn+2];
     7814      if (fabs(value) > zeroTolerance) {
     7815        array[numberNonZero] = value;
     7816        index[numberNonZero++] = column_[iColumn];
     7817        //index[numberNonZero++]=jColumn;
     7818      }
     7819      // jColumn = column_[iColumn+1];
     7820      value = 0.0;
     7821      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
     7822      for (j = start; j < end; j++) {
     7823        int iRow = row_[j];
     7824        value += pi[iRow] * element_[j];
     7825      }
     7826      //}
     7827    }
     7828    if (fabs(value) > zeroTolerance) {
     7829      array[numberNonZero] = value;
     7830      index[numberNonZero++] = column_[iColumn];
     7831      //index[numberNonZero++]=jColumn;
     7832    }
     7833  }
     7834  // do odd ones
     7835  for (int i=0;i<numberNonZero;i++) {
     7836    int iSequence = index[i];
     7837    double alpha;
     7838    double oldValue;
     7839    double value;
     7840    int iStatus = (statusArray[iSequence] & 3) - 1;
     7841    if (iStatus) {
     7842      double mult = multiplier[iStatus-1];
     7843      alpha = array[i] * mult;
     7844      if (alpha > 0.0) {
     7845        oldValue = reducedCost[iSequence] * mult;
     7846        value = oldValue - tentativeTheta * alpha;
     7847        if (value < dualT) {
     7848          value = oldValue - upperTheta * alpha;
     7849          if (value < dualT && alpha >= acceptablePivot) {
     7850            upperTheta = (oldValue - dualT) / alpha;
     7851          }
     7852          // add to list
     7853          spareArray[numberRemaining] = alpha * mult;
     7854          spareIndex[numberRemaining++] = iSequence;
     7855        }
     7856      }
     7857    }
     7858  }
     7859  int numberOld=numberNonZero;
     7860  int nMax=0;
     7861  for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
     7862    // C) Can do two at a time (if so put odd one into start_)
     7863    // D) can use switch
     7864    blockStruct * block = block_ + iBlock;
     7865    //int numberPrice = block->numberInBlock_;
     7866    int numberPrice = block->firstBasic_;
     7867    int nel = block->numberElements_;
     7868    const int * COIN_RESTRICT row = row_ + block->startElements_;
     7869    const double * COIN_RESTRICT element = element_ + block->startElements_;
     7870    const int * COIN_RESTRICT column = column_ + block->startIndices_;
     7871    int nBlock=numberPrice>>COIN_AVX2_SHIFT;
     7872    numberPrice -= roundDown(numberPrice);
     7873#if defined(HASWELL) || defined(SKYLAKE)
     7874    double * COIN_RESTRICT newValues = roundUpDouble((array+numberNonZero));
     7875#ifdef HASWELL
     7876    assert(COIN_AVX2==4);
     7877    for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7878      __m256d arrayX = _mm256_setzero_pd();
     7879      for (int j=0;j<nel;j++) {
     7880        __m128i rows = _mm_load_si128((const __m128i *)row); // was loadu
     7881        __m256d elements = _mm256_load_pd(element);
     7882        __m256d pis = _mm256_i32gather_pd(pi,rows,8);
     7883        arrayX = _mm256_fmadd_pd(pis,elements,arrayX);
     7884        row+=4;
     7885        element+=4;
     7886      }
     7887      _mm256_store_pd(newValues,arrayX);
     7888#else
     7889    assert(COIN_AVX2==8);
     7890    for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7891      __m512d arrayX = _mm512_setzero_pd();
     7892      for (int j=0;j<nel;j++) {
     7893        __m256i rows = _mm256_load_si256((const __m256i *)row); // was loadu
     7894        __m512d elements = _mm512_load_pd(element);
     7895        __m512d pis = _mm512_i32gather_pd(rows,pi,8);
     7896        arrayX = _mm512_fmadd_pd(pis,elements,arrayX);
     7897        row+=4;
     7898        element+=4;
     7899      }
     7900      _mm512_store_pd(newValues,arrayX);
     7901#endif
     7902      newValues += COIN_AVX2;
     7903      //row += (nel-1)*COIN_AVX2;
     7904      //element += (nel-1)*COIN_AVX2;
     7905      assert (row == row_ + block->startElements_ + nel*COIN_AVX2*(jBlock+1));
     7906    }
     7907#else
     7908    double * COIN_RESTRICT newValues = array+numberNonZero;
     7909    for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7910      for (int j=0; j < COIN_AVX2; j++) {
     7911        double value = 0.0;
     7912        for (int i=0; i < nel; i++) {
     7913          int iRow = row[i*COIN_AVX2];
     7914          value += pi[iRow] * element[i*COIN_AVX2];
     7915        }
     7916#if COIN_AVX2>1
     7917        row++;
     7918        element++;
     7919#else
     7920        row+=nel;
     7921        element+=nel;
     7922#endif
     7923        *newValues=value;
     7924        newValues++;
     7925      }
     7926#if COIN_AVX2>1
     7927      row += (nel-1)*COIN_AVX2;
     7928      element += (nel-1)*COIN_AVX2;
     7929      assert (row == row_ + block->startElements_ + nel*COIN_AVX2*(jBlock+1));
     7930#endif
     7931    }
     7932#endif
     7933    for (int j=0; j < numberPrice; j++) {
     7934      double value = 0.0;
     7935      for (int i=0; i < nel; i++) {
     7936        int iRow = row[i*COIN_AVX2];
     7937        value += pi[iRow] * element[i*COIN_AVX2];
     7938      }
     7939#if COIN_AVX2>1
     7940      row++;
     7941      element++;
     7942#else
     7943      row+=nel;
     7944      element+=nel;
     7945#endif
     7946      *newValues=value;
     7947      newValues++;
     7948    }
     7949#ifdef HASWELL
     7950    newValues = roundUpDouble((array+numberNonZero));
     7951#else
     7952    newValues = array+numberNonZero;
     7953#endif
     7954    int n = block->firstBasic_;
     7955    int nL = block->firstAtUpper_;
     7956    nMax = (newValues-array) + n;
     7957    for (int j=0;j<nL;j++) {
     7958      double value2=newValues[j];
     7959      if (fabs(value2) > zeroTolerance) {
     7960        int iSequence = column[j];
     7961        double alpha;
     7962        double oldValue;
     7963        double value;
     7964        alpha = value2 ;
     7965        if (alpha > 0.0) {
     7966          oldValue = reducedCost[iSequence];
     7967          value = oldValue - tentativeTheta * alpha;
     7968          if (value < dualT) {
     7969            value = oldValue - upperTheta * alpha;
     7970            if (value < dualT && alpha >= acceptablePivot) {
     7971              upperTheta = (oldValue - dualT) / alpha;
     7972            }
     7973            // add to list
     7974            spareArray[numberRemaining] = value2;
     7975            spareIndex[numberRemaining++] = iSequence;
     7976          }
     7977        }
     7978        array[numberNonZero] = value2;
     7979        index[numberNonZero++] = iSequence;
     7980      }
     7981    }
     7982    for (int j=nL;j<n;j++) {
     7983      double value2=newValues[j];
     7984      if (fabs(value2) > zeroTolerance) {
     7985        int iSequence = column[j];
     7986        double alpha;
     7987        double oldValue;
     7988        double value;
     7989        alpha = -value2;
     7990        if (alpha > 0.0) {
     7991          oldValue = -reducedCost[iSequence];
     7992          value = oldValue - tentativeTheta * alpha;
     7993          if (value < dualT) {
     7994            value = oldValue - upperTheta * alpha;
     7995            if (value < dualT && alpha >= acceptablePivot) {
     7996              upperTheta = (oldValue - dualT) / alpha;
     7997            }
     7998            // add to list
     7999            spareArray[numberRemaining] = value2;
     8000            spareIndex[numberRemaining++] = iSequence;
     8001          }
     8002        }
     8003        array[numberNonZero] = value2;
     8004        index[numberNonZero++] = iSequence;
     8005      }
     8006    }
     8007    numberOld=numberNonZero;
     8008  }
     8009  for (int j=numberNonZero;j<nMax;j++)
     8010    array[j]=0.0;
     8011  output->setNumElements(numberNonZero);
     8012  candidate->setNumElements(numberRemaining);
     8013  model->spareDoubleArray_[0] = upperTheta;
    76418014}
    76428015static void
     
    77668139  }
    77678140  info.numberAdded = bestSequence;
    7768   info.bestPossible = bestRatio;
    7769 }
    7770 #ifdef _MSC_VER
    7771 #include <intrin.h>
    7772 #else
    7773 #include <immintrin.h>
    7774 #endif
     8141}
    77758142
    77768143static void
     
    78178184#if COIN_AVX2 > 1
    78188185    int numberDo = roundDown(numberPrice);
     8186#if defined(HASWELL) || defined(SKYLAKE)
    78198187#ifdef HASWELL
    78208188    assert(COIN_AVX2==4);
     
    78268194        __m256d tempX = _mm256_setzero_pd();
    78278195        for (int j=0;j<nel;j++) {
    7828           __m128i rows = _mm_loadu_si128((const __m128i *)row);
     8196          __m128i rows = _mm_load_si128((const __m128i *)row); // was loadu
    78298197          __m256d elements = _mm256_load_pd(element);
    78308198          __m256d pis = _mm256_i32gather_pd(pi,rows,8);
     
    78408208        _mm256_store_pd(work2+i,arrayX);
    78418209      }
     8210#else
     8211    assert(COIN_AVX2==8);
     8212    __m512d zero = _mm512_setzero_pd();
     8213    for (int kColumn=0;kColumn<numberDo;kColumn+=COIN_AVX2_CHUNK) {
     8214      int endBlock=CoinMin(COIN_AVX2_CHUNK,numberDo-kColumn);
     8215      for(int i=0;i<endBlock;i+=COIN_AVX2) {
     8216        __m512d arrayX = _mm512_setzero_pd();
     8217        __m512d tempX = _mm512_setzero_pd();
     8218        for (int j=0;j<nel;j++) {
     8219          __m256i rows = _mm256_load_si256((const __m256i *)row); // was loadu
     8220          __m512d elements = _mm512_load_pd(element);
     8221          __m512d pis = _mm512_i32gather_pd(rows,pi,8);
     8222          __m512d piWeights = _mm512_i32gather_pd(rows,piWeight,8);
     8223          // should be better way using sub
     8224          arrayX = _mm512_fmadd_pd(pis,elements,arrayX);
     8225          tempX = _mm512_fmadd_pd(piWeights,elements,tempX);
     8226          row+=4;
     8227          element+=4;
     8228        }
     8229        arrayX = _mm512_sub_pd(zero,arrayX);
     8230        _mm512_store_pd(work+i,tempX);
     8231        _mm512_store_pd(work2+i,arrayX);
     8232      }
     8233#endif
    78428234      for (int i=0;i<endBlock;i++) {
    78438235        double value=work2[i];
     
    79308322  const blockStruct * block = blocks + iBlock;
    79318323  int first = 0;
    7932   int last = block->firstAtUpper_;
     8324  int last = block->firstAtLower_;
    79338325  const int * column = columnBlock + block->startIndices_;
    79348326  for (int j=first;j<last;j++) {
     
    79488340  }
    79498341  first = last;
    7950   last = block->firstAtLower_;
     8342  last = block->firstAtUpper_;
     8343  dualTolerance = - dualTolerance; //flip sign
    79518344  for (int j=first;j<last;j++) {
    79528345    int iColumn=*column;
    79538346    column++;
    79548347    double value = reducedCost[iColumn];
    7955     // at upper
    7956     if (value > dualTolerance) {
     8348    // at lower
     8349    if (value < dualTolerance) {
    79578350      value *= value;
    79588351      if (value>bestRatio*weights[iColumn]) {
     
    79698362    column++;
    79708363    double value = reducedCost[iColumn];
    7971     // at lower
    7972     if (value < dualTolerance) {
     8364    // at upper
     8365    if (value > dualTolerance) {
    79738366      value *= value;
    79748367      if (value>bestRatio*weights[iColumn]) {
     
    83928785#undef reference
    83938786#endif
    8394  
    8395  
     8787 #if 0
     8788#undef ABCSTATE_LITE
     8789#if ABOCA_LITE
     8790/* Meat of transposeTimes by column when not scaled and skipping
     8791   and doing part of dualColumn */
     8792static void
     8793dualColumn00(clpTempInfo & info)
     8794{
     8795  const int * COIN_RESTRICT which = info.which;
     8796  const double * COIN_RESTRICT work = info.work;
     8797  int * COIN_RESTRICT index = info.index;
     8798  double * COIN_RESTRICT spare = info.spare;
     8799  const unsigned char * COIN_RESTRICT status = info.status;
     8800  const double * COIN_RESTRICT reducedCost = info.reducedCost;
     8801  double upperTheta = info.upperTheta;
     8802  double acceptablePivot = info.acceptablePivot;
     8803  double dualTolerance = info.tolerance;
     8804  int numberToDo=info.numberToDo;
     8805  double tentativeTheta = 1.0e15;
     8806  int numberRemaining = 0;
     8807  double multiplier[] = { -1.0, 1.0};
     8808  double dualT = - dualTolerance;
     8809  for (int i = 0; i < numberToDo; i++) {
     8810    int iSequence = which[i];
     8811    int wanted = (status[iSequence] & 3) - 1;
     8812    if (wanted) {
     8813      double mult = multiplier[wanted-1];
     8814      double alpha = work[i] * mult;
     8815      if (alpha > 0.0) {
     8816        double oldValue = reducedCost[iSequence] * mult;
     8817        double value = oldValue - tentativeTheta * alpha;
     8818        if (value < dualT) {
     8819          value = oldValue - upperTheta * alpha;
     8820          if (value < dualT && alpha >= acceptablePivot) {
     8821            upperTheta = (oldValue - dualT) / alpha;
     8822          }
     8823          // add to list
     8824          spare[numberRemaining] = alpha * mult;
     8825          index[numberRemaining++] = iSequence;
     8826        }
     8827      }
     8828    }
     8829  }
     8830  info.numberRemaining = numberRemaining;
     8831  info.upperTheta = upperTheta;
     8832}
     8833#endif
     8834int
     8835ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
     8836                            const CoinIndexedVector * columnArray,
     8837                            CoinIndexedVector * spareArray,
     8838                            double acceptablePivot,
     8839                            double & upperReturn, double &bestReturn, double & badFree)
     8840{
     8841     // do first pass to get possibles
     8842     double * spare = spareArray->denseVector();
     8843     int * index = spareArray->getIndices();
     8844     const double * work;
     8845     int number;
     8846     const int * which;
     8847     const double * reducedCost;
     8848     // We can also see if infeasible or pivoting on free
     8849     double tentativeTheta = 1.0e15;
     8850     double upperTheta = 1.0e31;
     8851     double freePivot = acceptablePivot;
     8852     double bestPossible = 0.0;
     8853     int numberRemaining = 0;
     8854     int i;
     8855     badFree = 0.0;
     8856     if ((moreSpecialOptions_ & 8) != 0) {
     8857          // No free or super basic
     8858       // bestPossible will re recomputed if necessary
     8859       bestPossible=1.0;
     8860       //#define COIN_AVX7
     8861#ifndef COIN_AVX2
     8862          double multiplier[] = { -1.0, 1.0};
     8863#else
     8864          double multiplier[4]={0.0,0.0,-1.0,1.0};
     8865#endif
     8866          double dualT = - dualTolerance_;
     8867#if ABOCA_LITE==0
     8868          int nSections=2;
     8869#else
     8870          int numberThreads=abcState();
     8871          int nSections=numberThreads ? 1 : 2;
     8872#endif
     8873          for (int iSection = 0; iSection < nSections; iSection++) {
     8874
     8875               int addSequence;
     8876               unsigned char * statusArray;
     8877               if (!iSection) {
     8878                    work = rowArray->denseVector();
     8879                    number = rowArray->getNumElements();
     8880                    which = rowArray->getIndices();
     8881                    reducedCost = rowReducedCost_;
     8882                    addSequence = numberColumns_;
     8883                    statusArray = status_ + numberColumns_;
     8884               } else {
     8885                    work = columnArray->denseVector();
     8886                    number = columnArray->getNumElements();
     8887                    which = columnArray->getIndices();
     8888                    reducedCost = reducedCostWork_;
     8889                    addSequence = 0;
     8890                    statusArray = status_;
     8891               }
     8892#ifndef COIN_AVX2
     8893               for (i = 0; i < number; i++) {
     8894                    int iSequence = which[i];
     8895                    double alpha;
     8896                    double oldValue;
     8897                    double value;
     8898
     8899                    assert (getStatus(iSequence + addSequence) != isFree
     8900                            && getStatus(iSequence + addSequence) != superBasic);
     8901                    int iStatus = (statusArray[iSequence] & 3) - 1;
     8902                    if (iStatus) {
     8903                         double mult = multiplier[iStatus-1];
     8904                         alpha = work[i] * mult;
     8905                         if (alpha > 0.0) {
     8906                              oldValue = reducedCost[iSequence] * mult;
     8907                              value = oldValue - tentativeTheta * alpha;
     8908                              if (value < dualT) {
     8909                                   value = oldValue - upperTheta * alpha;
     8910                                   if (value < dualT && alpha >= acceptablePivot) {
     8911                                        upperTheta = (oldValue - dualT) / alpha;
     8912                                        //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     8913                                   }
     8914                                   // add to list
     8915                                   spare[numberRemaining] = alpha * mult;
     8916                                   index[numberRemaining++] = iSequence + addSequence;
     8917                              }
     8918                         }
     8919                    }
     8920               }
     8921               //
     8922#else
     8923               //#define COIN_AVX2 4 // temp
     8924#if COIN_AVX2 == 1
     8925#define COIN_AVX2_SHIFT 0
     8926#elif COIN_AVX2 == 2
     8927#define COIN_AVX2_SHIFT 1
     8928#elif COIN_AVX2 == 4
     8929#define COIN_AVX2_SHIFT 2
     8930#elif COIN_AVX2 == 8
     8931#define COIN_AVX2_SHIFT 3
     8932#else
     8933  error;
     8934#endif
     8935  //#define COIN_ALIGN 8*COIN_AVX2 // later
     8936  //#define COIN_ALIGN_DOUBLE COIN_AVX2
     8937#define CHECK_CHUNK 4
     8938               // round up
     8939               int * whichX = const_cast<int *>(which);
     8940               int nBlocks = (number+CHECK_CHUNK-1)/CHECK_CHUNK;
     8941               int n=nBlocks*CHECK_CHUNK+1;
     8942               for (int i=number;i<n;i++)
     8943                 whichX[i]=0; // alpha will be zero so not chosen
     8944               bool acceptableX[CHECK_CHUNK+1];
     8945               double oldValueX[CHECK_CHUNK+1];
     8946               double newValueX[CHECK_CHUNK+1];
     8947               double alphaX[CHECK_CHUNK+1];
     8948               newValueX[CHECK_CHUNK]=0.0;
     8949               #define USE_USE_AVX
     8950               //#define CHECK_H 1
     8951#ifdef USE_USE_AVX
     8952#define NEED_AVX
     8953#elif CHECK_H
     8954#define NEED_AVX
     8955#endif
     8956#ifdef NEED_AVX
     8957               double mult2[CHECK_CHUNK];
     8958               CoinInt64 acceptableY[CHECK_CHUNK];
     8959               CoinInt64 goodDj[CHECK_CHUNK];
     8960               double oldValueY[CHECK_CHUNK];
     8961               double alphaY[CHECK_CHUNK+1];
     8962               memset(acceptableY,0,sizeof(acceptableY));
     8963               memset(goodDj,0,sizeof(goodDj));
     8964               memset(oldValueY,0,sizeof(oldValueY));
     8965               memset(alphaY,0,sizeof(alphaY));
     8966               __m256d tentative2 = _mm256_set1_pd(-tentativeTheta);
     8967               __m256d dualT2 = _mm256_set1_pd(dualT);
     8968               __m256d acceptable2 = _mm256_set1_pd(acceptablePivot);
     8969#endif
     8970               for (int iBlock=0;iBlock<nBlocks;iBlock++) {
     8971                 bool store=false;
     8972                 double alpha=0.0;
     8973                 double oldValue=0.0;
     8974                 double newValue=0.0;
     8975                 double trueAlpha=0.0;
     8976                 int jSequence=0;
     8977#ifndef USE_USE_AVX
     8978                 for (int i = 0; i < CHECK_CHUNK+1; i++) {
     8979                   int iSequence = which[i];
     8980                   int iStatus = (statusArray[iSequence] & 3);
     8981                   double mult = multiplier[iStatus];
     8982                   double newAlpha = work[i]*mult;
     8983                   double oldDj = reducedCost[iSequence]*mult;
     8984                   newValue = (oldDj - tentativeTheta * newAlpha)-dualT;
     8985                   acceptableX[i]=newAlpha>= acceptablePivot;
     8986                   oldValueX[i]=oldDj;
     8987                   newValueX[i]=newValue;
     8988                   alphaX[i]=newAlpha;
     8989                 }
     8990#endif
     8991#ifdef NEED_AVX
     8992                 __m128i columns = _mm_load_si128((const __m128i *)which);
     8993                 // what do we get - this must be wrong
     8994                 // probably only 1 and 2 - can we be clever
     8995                 // fix
     8996                 //__m128i status; // = _mm256_i32gather_ps(statusArray,columns,1);
     8997                 //status.m128i_i32[0]=statusArray[columns.m128i_i32[0]];
     8998                 for (int i=0;i<CHECK_CHUNK;i++) {
     8999                   int iSequence = which[i];
     9000                   int iStatus = (statusArray[iSequence] & 3);
     9001                   mult2[i] = multiplier[iStatus];
     9002                 }
     9003                 //__m256d newAlpha2 = _mm256_i32gather_pd(multiplier,status,1); // mult here
     9004                 __m256d newAlpha2 = _mm256_load_pd(mult2);
     9005                 __m256d oldDj = _mm256_i32gather_pd(reducedCost,columns,8);
     9006                 oldDj=_mm256_mul_pd(oldDj,newAlpha2); // remember newAlpha==mult
     9007                 _mm256_store_pd(oldValueY,oldDj); // redo later
     9008                 __m256d work2 = _mm256_load_pd(work);
     9009                 newAlpha2=_mm256_mul_pd(newAlpha2,work2); // now really newAlpha
     9010                 //__m256d newValue2 = _mm256_fmadd_pd(tentative2,newAlpha2,oldDj);
     9011                 oldDj=_mm256_fmadd_pd(newAlpha2,tentative2,oldDj);
     9012                 __v4df bitsDj=_mm256_cmp_pd(oldDj,dualT2,_CMP_LT_OS);
     9013                 __v4df bitsAcceptable = _mm256_cmp_pd(newAlpha2,acceptable2,_CMP_GE_OS);
     9014                 _mm256_store_pd(reinterpret_cast<double *>(goodDj),bitsDj);
     9015                 _mm256_store_pd(reinterpret_cast<double *>(acceptableY),bitsAcceptable);
     9016                 _mm256_store_pd(alphaY,newAlpha2);
     9017#ifndef USE_USE_AVX
     9018#undef NDEBUG
     9019                 for (int i = 0; i < CHECK_CHUNK; i++) {
     9020                   assert(newValueX[i]>0.0==(goodDj[i]));
     9021                   //assert(acceptableX[i]==(acceptableY[i]));
     9022                   assert(oldValueX[i]==oldValueY[i]);
     9023                   assert(alphaX[i]==alphaY[i]);
     9024                 }
     9025                 for (int i = 0; i < CHECK_CHUNK; i++) {
     9026                   bool g1=newValueX[i]<0.0;
     9027                   bool g2=goodDj[i]!=0;
     9028                   if(g1!=g2)abort();
     9029                   //if(acceptableX[i]!=(acceptableY[i]))abort();
     9030                   if(fabs(oldValueX[i]-oldValueY[i])>1.0e-5+
     9031                      +(1.0e-10*fabs(oldValueX[i])))abort();
     9032                   if(alphaX[i]!=alphaY[i])abort();
     9033                 }
     9034#endif
     9035#endif
     9036                 for (int i = 0; i < CHECK_CHUNK+1; i++) {
     9037#ifndef USE_USE_AVX
     9038                   double newValue=newValueX[i];
     9039                   bool newStore = newValue < 0.0;
     9040                   if (store) {
     9041                     // add to list
     9042                     bool acceptable = acceptableX[i-1];
     9043                     spare[numberRemaining] = work[i-1];
     9044                     index[numberRemaining++] = which[i-1] + addSequence;
     9045                     double value = oldValueX[i-1] - upperTheta * alphaX[i-1];
     9046                     if (value < dualT && acceptable) {
     9047                       upperTheta = (oldValueX[i-1] - dualT) / alphaX[i-1];
     9048                     }
     9049                   }
     9050#else
     9051                   bool newStore = goodDj[i]!=0;
     9052                   if (store) {
     9053                     // add to list
     9054                     bool acceptable = acceptableY[i-1];
     9055                     spare[numberRemaining] = work[i-1];
     9056                     index[numberRemaining++] = which[i-1] + addSequence;
     9057                     double value = oldValueY[i-1] - upperTheta * alphaY[i-1];
     9058                     if (value < dualT && acceptable) {
     9059                       upperTheta = (oldValueY[i-1] - dualT) / alphaY[i-1];
     9060                     }
     9061                   }
     9062#endif
     9063                   store=newStore;
     9064                 }
     9065                 which += CHECK_CHUNK;
     9066                 work += CHECK_CHUNK;
     9067               }
     9068#endif
     9069          }
     9070#if ABOCA_LITE
     9071          if (numberThreads) {
     9072          work = columnArray->denseVector();
     9073          number = columnArray->getNumElements();
     9074          which = columnArray->getIndices();
     9075          reducedCost = reducedCostWork_;
     9076          unsigned char * statusArray = status_;
     9077         
     9078          clpTempInfo info[ABOCA_LITE];
     9079          int chunk = (number+numberThreads-1)/numberThreads;
     9080          int n=0;
     9081          int nR=numberRemaining;
     9082          for (int i=0;i<numberThreads;i++) {
     9083            info[i].which=const_cast<int *>(which+n);
     9084            info[i].work=const_cast<double *>(work+n);
     9085            info[i].numberToDo=CoinMin(chunk,number-n);
     9086            n += chunk;
     9087            info[i].index = index+nR;
     9088            info[i].spare = spare+nR;
     9089            nR += chunk;
     9090            info[i].reducedCost = const_cast<double *>(reducedCost);
     9091            info[i].upperTheta = upperTheta;
     9092            info[i].acceptablePivot = acceptablePivot;
     9093            info[i].status = statusArray;
     9094            info[i].tolerance=dualTolerance_;
     9095          }
     9096          for (int i=0;i<numberThreads;i++) {
     9097            cilk_spawn dualColumn00(info[i]);
     9098          }
     9099          cilk_sync;
     9100          moveAndZero(info,1,NULL);
     9101          for (int i=0;i<numberThreads;i++) {
     9102            numberRemaining += info[i].numberRemaining;
     9103            upperTheta = CoinMin(upperTheta,static_cast<double>(info[i].upperTheta));
     9104          }
     9105          }
     9106#endif
     9107     } else {
     9108          // some free or super basic
     9109          for (int iSection = 0; iSection < 2; iSection++) {
     9110
     9111               int addSequence;
     9112
     9113               if (!iSection) {
     9114                    work = rowArray->denseVector();
     9115                    number = rowArray->getNumElements();
     9116                    which = rowArray->getIndices();
     9117                    reducedCost = rowReducedCost_;
     9118                    addSequence = numberColumns_;
     9119               } else {
     9120                    work = columnArray->denseVector();
     9121                    number = columnArray->getNumElements();
     9122                    which = columnArray->getIndices();
     9123                    reducedCost = reducedCostWork_;
     9124                    addSequence = 0;
     9125               }
     9126
     9127               for (i = 0; i < number; i++) {
     9128                    int iSequence = which[i];
     9129                    double alpha;
     9130                    double oldValue;
     9131                    double value;
     9132                    bool keep;
     9133
     9134                    switch(getStatus(iSequence + addSequence)) {
     9135
     9136                    case basic:
     9137                    case ClpSimplex::isFixed:
     9138                         break;
     9139                    case isFree:
     9140                    case superBasic:
     9141                         alpha = work[i];
     9142                         bestPossible = CoinMax(bestPossible, fabs(alpha));
     9143                         oldValue = reducedCost[iSequence];
     9144                         // If free has to be very large - should come in via dualRow
     9145                         //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
     9146                         //break;
     9147                         if (oldValue > dualTolerance_) {
     9148                              keep = true;
     9149                         } else if (oldValue < -dualTolerance_) {
     9150                              keep = true;
     9151                         } else {
     9152                              if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
     9153                                   keep = true;
     9154                              } else {
     9155                                   keep = false;
     9156                                   badFree = CoinMax(badFree, fabs(alpha));
     9157                              }
     9158                         }
     9159                         if (keep) {
     9160                              // free - choose largest
     9161                              if (fabs(alpha) > freePivot) {
     9162                                   freePivot = fabs(alpha);
     9163                                   sequenceIn_ = iSequence + addSequence;
     9164                                   theta_ = oldValue / alpha;
     9165                                   alpha_ = alpha;
     9166                              }
     9167                              // give fake bounds if possible
     9168                              int jSequence=iSequence+addSequence;
     9169                              if (2.0*fabs(solution_[jSequence])<
     9170                                  dualBound_) {
     9171                                FakeBound bound = getFakeBound(jSequence);
     9172                                assert (bound == ClpSimplexDual::noFake);
     9173                                setFakeBound(jSequence,ClpSimplexDual::bothFake);
     9174                                numberFake_++;
     9175                                value = oldValue - tentativeTheta * alpha;
     9176                                if (value > dualTolerance_) {
     9177                                  // pretend coming in from upper bound
     9178                                  upper_[jSequence] = solution_[jSequence];
     9179                                  lower_[jSequence] = upper_[jSequence] - dualBound_;
     9180                                  setColumnStatus(jSequence,ClpSimplex::atUpperBound);
     9181                                } else {
     9182                                  // pretend coming in from lower bound
     9183                                  lower_[jSequence] = solution_[jSequence];
     9184                                  upper_[jSequence] = lower_[jSequence] + dualBound_;
     9185                                  setColumnStatus(jSequence,ClpSimplex::atLowerBound);
     9186                                }
     9187                              }
     9188                         }
     9189                         break;
     9190                    case atUpperBound:
     9191                         alpha = work[i];
     9192                         oldValue = reducedCost[iSequence];
     9193                         value = oldValue - tentativeTheta * alpha;
     9194                         //assert (oldValue<=dualTolerance_*1.0001);
     9195                         if (value > dualTolerance_) {
     9196                              bestPossible = CoinMax(bestPossible, -alpha);
     9197                              value = oldValue - upperTheta * alpha;
     9198                              if (value > dualTolerance_ && -alpha >= acceptablePivot) {
     9199                                   upperTheta = (oldValue - dualTolerance_) / alpha;
     9200                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     9201                              }
     9202                              // add to list
     9203                              spare[numberRemaining] = alpha;
     9204                              index[numberRemaining++] = iSequence + addSequence;
     9205                         }
     9206                         break;
     9207                    case atLowerBound:
     9208                         alpha = work[i];
     9209                         oldValue = reducedCost[iSequence];
     9210                         value = oldValue - tentativeTheta * alpha;
     9211                         //assert (oldValue>=-dualTolerance_*1.0001);
     9212                         if (value < -dualTolerance_) {
     9213                              bestPossible = CoinMax(bestPossible, alpha);
     9214                              value = oldValue - upperTheta * alpha;
     9215                              if (value < -dualTolerance_ && alpha >= acceptablePivot) {
     9216                                   upperTheta = (oldValue + dualTolerance_) / alpha;
     9217                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     9218                              }
     9219                              // add to list
     9220                              spare[numberRemaining] = alpha;
     9221                              index[numberRemaining++] = iSequence + addSequence;
     9222                         }
     9223                         break;
     9224                    }
     9225               }
     9226          }
     9227     }
     9228     upperReturn = upperTheta;
     9229     bestReturn = bestPossible;
     9230     return numberRemaining;
     9231}
     9232#endif
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r2235 r2259  
    410410                                      const double * COIN_RESTRICT reducedCost,
    411411                                      double & upperTheta,
    412                                       double & bestPossible,
    413412                                      double acceptablePivot,
    414413                                      double dualTolerance,
     
    488487     int * indexTemp;
    489488     int * numberInPtr;
    490      double * bestPossiblePtr;
     489  //double * bestPossiblePtr;
    491490     double * upperThetaPtr;
    492491     int * posFreePtr;
     
    573572  /* order is -
    574573     free or superbasic
     574     at lower
    575575     at upper
    576      at lower
    577576     fixed or basic */
     577  int firstAtLower_;
    578578  int firstAtUpper_;
    579   int firstAtLower_;
    580579  int firstBasic_; // or fixed
    581580  int numberElements_; // number elements per column
     
    593592                         const double * pi,
    594593                         CoinIndexedVector * output) const;
     594     /// This version does dualColumn0
    595595     /// Updates two arrays for steepest
     596     void transposeTimes(const ClpSimplex * model,
     597                         const double * pi,
     598                         CoinIndexedVector * output,
     599                         CoinIndexedVector * candidate,
     600                         const CoinIndexedVector * rowArray) const;
    596601     void transposeTimes2(const ClpSimplex * model,
    597602                          const double * pi, CoinIndexedVector * dj1,
     
    631636     /// Part of above
    632637     void swapOne(int iBlock, int kA, int kB);
     638     /** Debug - check blocks */
     639     void checkBlocks(const ClpSimplex * model);
    633640     /**
    634641        type - 1 redo infeasible, 2 choose sequenceIn, 3 both
     
    676683     /// Blocks (ordinary start at 0 and go to first block)
    677684     blockStruct * block_;
     685     /// If active
     686     int ifActive_;
    678687     //@}
    679688};
  • trunk/Clp/src/ClpPresolve.cpp

    r2235 r2259  
    22592259          if ((presolveActions_ & 0x80000000) != 0)
    22602260               prob.statistics();
     2261          if (doTransfer())
     2262              transferCosts(&prob);
    22612263          // make sure row solution correct
    22622264          {
  • trunk/Clp/src/ClpPresolve.hpp

    r2146 r2259  
    160160          else presolveActions_ &= ~32768;
    161161     }
     162     /// Whether we want to do transfer part of presolve
     163     inline bool doTransfer() const {
     164          return (presolveActions_ & 65536) != 0;
     165     }
     166     inline void setDoTransfer(bool doTransfer) {
     167          if (doTransfer) presolveActions_  |= 65536;
     168          else presolveActions_ &= ~65536;
     169     }
    162170     /// Whether we want to do singleton column part of presolve
    163171     inline bool doSingletonColumn() const {
     
    203211     /// Set whole group
    204212     inline int presolveActions() const {
    205           return presolveActions_ & 0xffff;
     213          return presolveActions_ & 0xffffff;
    206214     }
    207215     inline void setPresolveActions(int action) {
    208           presolveActions_  = (presolveActions_ & 0xffff0000) | (action & 0xffff);
     216          presolveActions_  = (presolveActions_ & 0xff000000) | (action & 0xffffff);
    209217     }
    210218     /// Substitution level
  • trunk/Clp/src/ClpSimplex.cpp

    r2235 r2259  
    5353     moreSpecialOptions_(2),
    5454     baseIteration_(0),
     55     vectorMode_(0),
    5556     primalToleranceToGetOptimal_(-1.0),
    5657     largeValue_(1.0e15),
     
    168169     moreSpecialOptions_(2),
    169170     baseIteration_(0),
     171     vectorMode_(0),
    170172     primalToleranceToGetOptimal_(-1.0),
    171173     largeValue_(1.0e15),
     
    320322     moreSpecialOptions_(2),
    321323     baseIteration_(0),
     324     vectorMode_(0),
    322325     primalToleranceToGetOptimal_(-1.0),
    323326     largeValue_(1.0e15),
     
    23312334     moreSpecialOptions_(2),
    23322335     baseIteration_(0),
     2336     vectorMode_(0),
    23332337     primalToleranceToGetOptimal_(-1.0),
    23342338     largeValue_(1.0e15),
     
    24392443     moreSpecialOptions_(2),
    24402444     baseIteration_(0),
     2445     vectorMode_(0),
    24412446     primalToleranceToGetOptimal_(-1.0),
    24422447     largeValue_(1.0e15),
     
    26422647     bestObjectiveValue_ = rhs.bestObjectiveValue_;
    26432648     baseIteration_ = rhs.baseIteration_;
     2649     vectorMode_ = rhs.vectorMode_;
    26442650     primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
    26452651     largeValue_ = rhs.largeValue_;
     
    62516257          barrier.setCholesky(cholesky);
    62526258     }
    6253 #elif defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD)
     6259#elif 0 //defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD)
    62546260     if (!doKKT) {
    62556261          ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
     
    92189224     otherModel.moreSpecialOptions_ = moreSpecialOptions_;
    92199225     otherModel.baseIteration_ = baseIteration_;
     9226     otherModel.vectorMode_ = vectorMode_;
    92209227     otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
    92219228     otherModel.largestPrimalError_ = largestPrimalError_;
  • trunk/Clp/src/ClpSimplex.hpp

    r2235 r2259  
    1818#include "ClpSolve.hpp"
    1919#include "ClpConfig.h"
     20#include "CoinIndexedVector.hpp"
    2021class ClpDualRowPivot;
    2122class ClpPrimalColumnPivot;
     
    12881289          return moreSpecialOptions_;
    12891290     }
     1291     /// Get vector mode
     1292     inline int vectorMode() const {
     1293          return vectorMode_;
     1294     }
    12901295     /** Set more special options
    12911296         1 bit - if presolve says infeasible in ClpSolve return
     
    13121317         4194304 bit - tolerances have been changed by code
    13131318         8388608 bit - tolerances are dynamic (at first)
    1314          16777216 bit - try vector matrix
    13151319     */
    13161320     inline void setMoreSpecialOptions(int value) {
    13171321          moreSpecialOptions_ = value;
     1322     }
     1323     /// Set vector mode
     1324     inline void setVectorMode(int value) {
     1325          vectorMode_ = value;
    13181326     }
    13191327     //@}
     
    15831591     /// Iteration when we entered dual or primal
    15841592     int baseIteration_;
     1593     /// Vector mode - try and use vector instructions
     1594     int vectorMode_;
    15851595     /// Primal tolerance needed to make dual feasible (<largeTolerance)
    15861596     double primalToleranceToGetOptimal_;
  • trunk/Clp/src/ClpSimplexDual.cpp

    r2235 r2259  
    631631          returnCode = 1; // to skip gutsOfDual
    632632          problemStatus_ = 0;
     633     } else if (maximumIterations()==0) {
     634          returnCode = 1; // to skip gutsOfDual
     635          problemStatus_ = 3;
    633636     }
    634637
     
    33813384  double acceptablePivot = info.acceptablePivot;
    33823385  double dualTolerance = info.tolerance;
    3383   double bestPossible = info.bestPossible;
    33843386  int numberToDo=info.numberToDo;
    33853387  double tentativeTheta = 1.0e15;
     
    33973399        double value = oldValue - tentativeTheta * alpha;
    33983400        if (value < dualT) {
    3399           bestPossible = CoinMax(bestPossible, alpha);
    34003401          value = oldValue - upperTheta * alpha;
    34013402          if (value < dualT && alpha >= acceptablePivot) {
     
    34113412  info.numberRemaining = numberRemaining;
    34123413  info.upperTheta = upperTheta;
    3413   info.bestPossible = bestPossible;
     3414}
     3415static void
     3416dualColumn000(int numberThreads,clpTempInfo * info)
     3417{
     3418  for (int i=0;i<numberThreads;i++) {
     3419    cilk_spawn dualColumn00(info[i]);
     3420  }
     3421  cilk_sync;
    34143422}
    34153423void moveAndZero(clpTempInfo * info,int type,void * extra)
     
    34833491}
    34843492#endif
     3493#ifdef _MSC_VER
     3494#include <intrin.h>
     3495#else
     3496#include <immintrin.h>
     3497//#include <fmaintrin.h>
     3498#endif
    34853499int
    34863500ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
     
    34883502                            CoinIndexedVector * spareArray,
    34893503                            double acceptablePivot,
    3490                             double & upperReturn, double &bestReturn, double & badFree)
     3504                            double & upperReturn, double & badFree)
    34913505{
    34923506     // do first pass to get possibles
     
    35013515     double upperTheta = 1.0e31;
    35023516     double freePivot = acceptablePivot;
    3503      double bestPossible = 0.0;
    35043517     int numberRemaining = 0;
    35053518     int i;
     
    35073520     if ((moreSpecialOptions_ & 8) != 0) {
    35083521          // No free or super basic
     3522       // bestPossible will re recomputed if necessary
     3523#ifndef COIN_AVX2
    35093524          double multiplier[] = { -1.0, 1.0};
     3525#else
     3526          double multiplier[4]={0.0,0.0,-1.0,1.0};
     3527#endif
    35103528          double dualT = - dualTolerance_;
    35113529#if ABOCA_LITE==0
     
    35343552                    statusArray = status_;
    35353553               }
    3536 
     3554#ifndef COIN_AVX2
    35373555               for (i = 0; i < number; i++) {
    35383556                    int iSequence = which[i];
     
    35513569                              value = oldValue - tentativeTheta * alpha;
    35523570                              if (value < dualT) {
    3553                                    bestPossible = CoinMax(bestPossible, alpha);
    35543571                                   value = oldValue - upperTheta * alpha;
    35553572                                   if (value < dualT && alpha >= acceptablePivot) {
     
    35643581                    }
    35653582               }
     3583               //
     3584#else
     3585               //#define COIN_AVX2 4 // temp
     3586#if COIN_AVX2 == 1
     3587#define COIN_AVX2_SHIFT 0
     3588#elif COIN_AVX2 == 2
     3589#define COIN_AVX2_SHIFT 1
     3590#elif COIN_AVX2 == 4
     3591#define COIN_AVX2_SHIFT 2
     3592#elif COIN_AVX2 == 8
     3593#define COIN_AVX2_SHIFT 3
     3594#else
     3595  error;
     3596#endif
     3597  //#define COIN_ALIGN 8*COIN_AVX2 // later
     3598  //#define COIN_ALIGN_DOUBLE COIN_AVX2
     3599#define CHECK_CHUNK 4
     3600               // round up
     3601               int * whichX = const_cast<int *>(which);
     3602               int nBlocks = (number+CHECK_CHUNK-1)/CHECK_CHUNK;
     3603               int n=nBlocks*CHECK_CHUNK+1;
     3604               for (int i=number;i<n;i++)
     3605                 whichX[i]=0; // alpha will be zero so not chosen
     3606               bool acceptableX[CHECK_CHUNK+1];
     3607               double oldValueX[CHECK_CHUNK+1];
     3608               double newValueX[CHECK_CHUNK+1];
     3609               double alphaX[CHECK_CHUNK+1];
     3610               newValueX[CHECK_CHUNK]=0.0;
     3611               #define USE_USE_AVX
     3612               //#define CHECK_H 1
     3613#ifdef USE_USE_AVX
     3614#define NEED_AVX
     3615#elif CHECK_H
     3616#define NEED_AVX
     3617#endif
     3618#ifdef NEED_AVX
     3619               double mult2[CHECK_CHUNK] __attribute__((aligned(64)));
     3620               CoinInt64 acceptableY[CHECK_CHUNK] __attribute__((aligned(64)));
     3621               CoinInt64 goodDj[CHECK_CHUNK] __attribute__((aligned(64)));
     3622               double oldValueY[CHECK_CHUNK] __attribute__((aligned(64)));
     3623               double alphaY[CHECK_CHUNK+1] __attribute__((aligned(64)));
     3624               memset(acceptableY,0,sizeof(acceptableY));
     3625               memset(goodDj,0,sizeof(goodDj));
     3626               memset(oldValueY,0,sizeof(oldValueY));
     3627               memset(alphaY,0,sizeof(alphaY));
     3628               __m256d tentative2 = _mm256_set1_pd(-tentativeTheta);
     3629               __m256d dualT2 = _mm256_set1_pd(dualT);
     3630               __m256d acceptable2 = _mm256_set1_pd(acceptablePivot);
     3631#endif
     3632               for (int iBlock=0;iBlock<nBlocks;iBlock++) {
     3633                 bool store=false;
     3634                 double alpha=0.0;
     3635                 double oldValue=0.0;
     3636                 double newValue=0.0;
     3637                 double trueAlpha=0.0;
     3638                 int jSequence=0;
     3639#ifndef USE_USE_AVX
     3640                 for (int i = 0; i < CHECK_CHUNK+1; i++) {
     3641                   int iSequence = which[i];
     3642                   int iStatus = (statusArray[iSequence] & 3);
     3643                   double mult = multiplier[iStatus];
     3644                   double newAlpha = work[i]*mult;
     3645                   double oldDj = reducedCost[iSequence]*mult;
     3646                   newValue = (oldDj - tentativeTheta * newAlpha)-dualT;
     3647                   acceptableX[i]=newAlpha>= acceptablePivot;
     3648                   oldValueX[i]=oldDj;
     3649                   newValueX[i]=newValue;
     3650                   alphaX[i]=newAlpha;
     3651                 }
     3652#endif
     3653#ifdef NEED_AVX
     3654                 __m128i columns = _mm_load_si128((const __m128i *)which);
     3655                 // what do we get - this must be wrong
     3656                 // probably only 1 and 2 - can we be clever
     3657                 // fix
     3658                 //__m128i status; // = _mm256_i32gather_ps(statusArray,columns,1);
     3659                 //status.m128i_i32[0]=statusArray[columns.m128i_i32[0]];
     3660                 for (int i=0;i<CHECK_CHUNK;i++) {
     3661                   int iSequence = which[i];
     3662                   int iStatus = (statusArray[iSequence] & 3);
     3663                   mult2[i] = multiplier[iStatus];
     3664                 }
     3665                 //__m256d newAlpha2 = _mm256_i32gather_pd(multiplier,status,1); // mult here
     3666                 __m256d newAlpha2 = _mm256_load_pd(mult2);
     3667                 __m256d oldDj = _mm256_i32gather_pd(reducedCost,columns,8);
     3668                 oldDj=_mm256_mul_pd(oldDj,newAlpha2); // remember newAlpha==mult
     3669                 _mm256_store_pd(oldValueY,oldDj); // redo later
     3670                 __m256d work2 = _mm256_load_pd(work);
     3671                 newAlpha2=_mm256_mul_pd(newAlpha2,work2); // now really newAlpha
     3672                 //__m256d newValue2 = _mm256_fmadd_pd(tentative2,newAlpha2,oldDj);
     3673                 oldDj=_mm256_fmadd_pd(newAlpha2,tentative2,oldDj);
     3674                 __v4df bitsDj=_mm256_cmp_pd(oldDj,dualT2,_CMP_LT_OS);
     3675                 __v4df bitsAcceptable = _mm256_cmp_pd(newAlpha2,acceptable2,_CMP_GE_OS);
     3676                 _mm256_store_pd(reinterpret_cast<double *>(goodDj),bitsDj);
     3677                 _mm256_store_pd(reinterpret_cast<double *>(acceptableY),bitsAcceptable);
     3678                 _mm256_store_pd(alphaY,newAlpha2);
     3679#ifndef USE_USE_AVX
     3680#undef NDEBUG
     3681                 for (int i = 0; i < CHECK_CHUNK; i++) {
     3682                   assert(newValueX[i]>0.0==(goodDj[i]));
     3683                   //assert(acceptableX[i]==(acceptableY[i]));
     3684                   assert(oldValueX[i]==oldValueY[i]);
     3685                   assert(alphaX[i]==alphaY[i]);
     3686                 }
     3687                 for (int i = 0; i < CHECK_CHUNK; i++) {
     3688                   bool g1=newValueX[i]<0.0;
     3689                   bool g2=goodDj[i]!=0;
     3690                   if(g1!=g2)abort();
     3691                   //if(acceptableX[i]!=(acceptableY[i]))abort();
     3692                   if(fabs(oldValueX[i]-oldValueY[i])>1.0e-5+
     3693                      +(1.0e-10*fabs(oldValueX[i])))abort();
     3694                   if(alphaX[i]!=alphaY[i])abort();
     3695                 }
     3696#endif
     3697#endif
     3698                 for (int i = 0; i < CHECK_CHUNK+1; i++) {
     3699#ifndef USE_USE_AVX
     3700                   double newValue=newValueX[i];
     3701                   bool newStore = newValue < 0.0;
     3702                   if (store) {
     3703                     // add to list
     3704                     bool acceptable = acceptableX[i-1];
     3705                     spare[numberRemaining] = work[i-1];
     3706                     index[numberRemaining++] = which[i-1] + addSequence;
     3707                     double value = oldValueX[i-1] - upperTheta * alphaX[i-1];
     3708                     if (value < dualT && acceptable) {
     3709                       upperTheta = (oldValueX[i-1] - dualT) / alphaX[i-1];
     3710                     }
     3711                   }
     3712#else
     3713                   bool newStore = goodDj[i]!=0;
     3714                   if (store) {
     3715                     // add to list
     3716                     bool acceptable = acceptableY[i-1];
     3717                     spare[numberRemaining] = work[i-1];
     3718                     index[numberRemaining++] = which[i-1] + addSequence;
     3719                     double value = oldValueY[i-1] - upperTheta * alphaY[i-1];
     3720                     if (value < dualT && acceptable) {
     3721                       upperTheta = (oldValueY[i-1] - dualT) / alphaY[i-1];
     3722                     }
     3723                   }
     3724#endif
     3725                   store=newStore;
     3726                 }
     3727                 which += CHECK_CHUNK;
     3728                 work += CHECK_CHUNK;
     3729               }
     3730#endif
    35663731          }
    35673732#if ABOCA_LITE
     
    35873752            info[i].reducedCost = const_cast<double *>(reducedCost);
    35883753            info[i].upperTheta = upperTheta;
    3589             info[i].bestPossible = bestPossible;
    35903754            info[i].acceptablePivot = acceptablePivot;
    35913755            info[i].status = statusArray;
    35923756            info[i].tolerance=dualTolerance_;
    35933757          }
    3594           for (int i=0;i<numberThreads;i++) {
    3595             cilk_spawn dualColumn00(info[i]);
    3596           }
    3597           cilk_sync;
     3758          // for gcc - get cilk out of function to stop avx2 error
     3759          dualColumn000(numberThreads,info);
    35983760          moveAndZero(info,1,NULL);
    35993761          for (int i=0;i<numberThreads;i++) {
    36003762            numberRemaining += info[i].numberRemaining;
    3601             bestPossible = CoinMax(bestPossible,static_cast<double>(info[i].bestPossible));
    36023763            upperTheta = CoinMin(upperTheta,static_cast<double>(info[i].upperTheta));
    36033764          }
     
    36393800                    case superBasic:
    36403801                         alpha = work[i];
    3641                          bestPossible = CoinMax(bestPossible, fabs(alpha));
    36423802                         oldValue = reducedCost[iSequence];
    36433803                         // If free has to be very large - should come in via dualRow
     
    36933853                         //assert (oldValue<=dualTolerance_*1.0001);
    36943854                         if (value > dualTolerance_) {
    3695                               bestPossible = CoinMax(bestPossible, -alpha);
    36963855                              value = oldValue - upperTheta * alpha;
    36973856                              if (value > dualTolerance_ && -alpha >= acceptablePivot) {
     
    37103869                         //assert (oldValue>=-dualTolerance_*1.0001);
    37113870                         if (value < -dualTolerance_) {
    3712                               bestPossible = CoinMax(bestPossible, alpha);
    37133871                              value = oldValue - upperTheta * alpha;
    37143872                              if (value < -dualTolerance_ && alpha >= acceptablePivot) {
     
    37263884     }
    37273885     upperReturn = upperTheta;
    3728      bestReturn = bestPossible;
    37293886     return numberRemaining;
    37303887}
     
    38143971     // do first pass to get possibles
    38153972     upperTheta = 1.0e31;
    3816      double bestPossible = 0.0;
     3973     double bestPossible = 1.0;
    38173974     double badFree = 0.0;
    38183975     alpha_ = 0.0;
    38193976     if (spareIntArray_[0] >= 0) {
    38203977          numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
    3821                                         acceptablePivot, upperTheta, bestPossible, badFree);
     3978                                        acceptablePivot, upperTheta, badFree);
    38223979     } else {
    38233980          // already done
     
    38253982          spareArray->setNumElements(0);
    38263983          upperTheta = spareDoubleArray_[0];
    3827           bestPossible = spareDoubleArray_[1];
    38283984          if (spareIntArray_[0] == -1) {
    38293985               theta_ = spareDoubleArray_[2];
     
    38323988          } else {
    38333989#if 0
     3990#undef NDEBUG
    38343991               int n = numberRemaining;
    38353992               double u = upperTheta;
    3836                double b = bestPossible;
    38373993               upperTheta = 1.0e31;
    3838                bestPossible = 0.0;
    3839                numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
    3840                                              acceptablePivot, upperTheta, bestPossible, badFree);
     3994               CoinIndexedVector temp(4000);
     3995               numberRemaining = dualColumn0(rowArray, columnArray, &temp,
     3996                                             acceptablePivot, upperTheta, badFree);
    38413997               assert (n == numberRemaining);
    3842                assert (fabs(b - bestPossible) < 1.0e-7);
     3998               double * spare = spareArray->denseVector();
     3999               int * index = spareArray->getIndices();
     4000               double * spareX = temp.denseVector();
     4001               int * indexX = temp.getIndices();
     4002               CoinSort_2(spare,spare+n,index);
     4003               CoinSort_2(spareX,spareX+n,indexX);
     4004               for (int i=0;i<n;i++) {
     4005                 assert (index[i]==indexX[i]);
     4006                 assert (fabs(spare[i]-spareX[i])<1.0e-6);
     4007               }
    38434008               assert (fabs(u - upperTheta) < 1.0e-7);
    38444009#endif
     
    44204585               lowerIn_ = valueIn_;
    44214586          }
     4587          if (fabs(alpha_)<1.0e-6) {
     4588            // need bestPossible
     4589            const double * work;
     4590            int number;
     4591            const int * which;
     4592            const double * reducedCost;
     4593            double tentativeTheta = 1.0e15;
     4594            double upperTheta = 1.0e31;
     4595            bestPossible = 0.0;
     4596            double multiplier[] = { -1.0, 1.0};
     4597            double dualT = - dualTolerance_;
     4598            int nSections=2;
     4599            int addSequence;
     4600            for (int iSection = 0; iSection < nSections; iSection++) {
     4601              if (!iSection) {
     4602                work = rowArray->denseVector();
     4603                number = rowArray->getNumElements();
     4604                which = rowArray->getIndices();
     4605                reducedCost = rowReducedCost_;
     4606                addSequence = numberColumns_;
     4607              } else {
     4608                work = columnArray->denseVector();
     4609                number = columnArray->getNumElements();
     4610                which = columnArray->getIndices();
     4611                reducedCost = reducedCostWork_;
     4612                addSequence = 0;
     4613              }
     4614              for (i = 0; i < number; i++) {
     4615                int iSequence = which[i];
     4616                double alpha;
     4617                double oldValue;
     4618                double value;
     4619                double mult=1.0;
     4620                switch(getStatus(iSequence + addSequence)) {
     4621                 
     4622                case basic:
     4623                case ClpSimplex::isFixed:
     4624                  break;
     4625                case isFree:
     4626                case superBasic:
     4627                  alpha = work[i];
     4628                  bestPossible = CoinMax(bestPossible, fabs(alpha));
     4629                  break;
     4630                case atUpperBound:
     4631                  mult = -1.0;
     4632                case atLowerBound:
     4633                  alpha = work[i] * mult;
     4634                  if (alpha > 0.0) {
     4635                    oldValue = reducedCost[iSequence] * mult;
     4636                    value = oldValue - tentativeTheta * alpha;
     4637                    if (value < dualT) {
     4638                      bestPossible = CoinMax(bestPossible, alpha);
     4639                    }
     4640                  }
     4641                  break;
     4642                }
     4643              }
     4644            }
     4645          } else {
     4646            bestPossible=fabs(alpha_);
     4647          }
    44224648     } else {
    44234649          // no pivot
     
    58836109               case atUpperBound:
    58846110                    // to lower bound
    5885                     setStatus(iSequence + addSequence, atLowerBound);
     6111                    setStatus(iSequence+addSequence, atLowerBound);
    58866112                    solution[iSequence] = lower[iSequence];
     6113                    // correct in vector copy
     6114                    iSequence += addSequence;
     6115                    matrix_->correctSequence(this,iSequence,iSequence);
    58876116                    break;
    58886117               case atLowerBound:
    58896118                    // to upper bound
    5890                     setStatus(iSequence + addSequence, atUpperBound);
     6119                    setStatus(iSequence+addSequence, atUpperBound);
    58916120                    solution[iSequence] = upper[iSequence];
     6121                    // correct in vector copy
     6122                    iSequence += addSequence;
     6123                    matrix_->correctSequence(this,iSequence,iSequence);
    58926124                    break;
    58936125               }
     
    59026134{
    59036135     if (getFakeBound(iSequence) != noFake) {
    5904           numberFake_--;;
     6136          numberFake_--;
    59056137          setFakeBound(iSequence, noFake);
    59066138          if (iSequence >= numberColumns_) {
  • trunk/Clp/src/ClpSimplexDual.hpp

    r1761 r2259  
    201201                     CoinIndexedVector * spareArray,
    202202                     double acceptablePivot,
    203                      double & upperReturn, double &bestReturn, double & badFree);
     203                     double & upperReturn, double & badFree);
    204204     /**
    205205         Row array has row part of pivot row
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r2254 r2259  
    824824     int dummy; // for use in generalExpanded
    825825     int saveFirstFree = firstFree_;
     826     double saveOriginalWeight = infeasibilityCost_;
    826827     // number of pivots done
    827828     int numberPivots = factorization_->pivots();
     
    13721373                    nonLinearCost_->checkInfeasibilities(primalTolerance_);
    13731374               } else {
     1375#ifndef MAX_INFEASIBILITY_COST
     1376#define MAX_INFEASIBILITY_COST 1.0e18
     1377#endif
    13741378                    infeasibilityCost_ = 1.0e30;
    13751379                    gutsOfSolution(NULL, NULL, ifValuesPass != 0 && firstFree_>=0);
    1376                     infeasibilityCost_ = saveWeight;
    1377                     if ((infeasibilityCost_ >= 1.0e18 ||
     1380                    infeasibilityCost_ = CoinMax(saveWeight,saveOriginalWeight);
     1381                    if ((infeasibilityCost_ >= MAX_INFEASIBILITY_COST ||
    13781382                              numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
    13791383                         goToDual = unPerturb(); // stop any further perturbation
     
    13891393                    }
    13901394                    if (!goToDual) {
    1391                          if (infeasibilityCost_ >= 1.0e20 ||
     1395                         if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST ||
    13921396                                   numberDualInfeasibilities_ == 0) {
    13931397                              // we are infeasible - use as ray
     
    14121416                                   nonLinearCost_->numberInfeasibilities();
    14131417                         }
    1414                          if (infeasibilityCost_ < 1.0e20) {
     1418                         if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
    14151419                              infeasibilityCost_ *= 5.0;
    14161420                              // reset looping criterion
     
    15291533          if (problemStatus_ == -5) {
    15301534               if (nonLinearCost_->numberInfeasibilities()) {
    1531                     if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) {
     1535                    if (infeasibilityCost_ >= MAX_INFEASIBILITY_COST && perturbation_ == 101) {
    15321536                         // back off weight
    15331537                         infeasibilityCost_ = 1.0e13;
     
    15371541                    }
    15381542                    //we need infeasiblity cost changed
    1539                     if (infeasibilityCost_ < 1.0e20) {
     1543                    if (infeasibilityCost_ < MAX_INFEASIBILITY_COST) {
    15401544                         infeasibilityCost_ *= 5.0;
    15411545                         // reset looping criterion
  • trunk/Clp/src/ClpSolve.cpp

    r2241 r2259  
    913913     bool presolveToFile = (presolveOptions & 0x40000000) != 0;
    914914     presolveOptions &= ~0x40000000;
    915      if ((presolveOptions & 0xffff) != 0)
     915     if ((presolveOptions & 0xffffff) != 0)
    916916          pinfo->setPresolveActions(presolveOptions);
    917917     // switch off singletons to slacks
     
    24502450
    24512451          // See if we have costed slacks
    2452           int * negSlack = new int[numberRows];
    2453           int * posSlack = new int[numberRows];
     2452          int * negSlack = new int[numberRows+numberColumns];
     2453          int * posSlack = new int[numberRows+numberColumns];
     2454          int * nextNegSlack = negSlack+numberRows;
     2455          int * nextPosSlack = posSlack+numberRows;
    24542456          int iRow;
    2455           for (iRow = 0; iRow < numberRows; iRow++) {
     2457          for (iRow = 0; iRow < numberRows+numberColumns; iRow++) {
    24562458               negSlack[iRow] = -1;
    24572459               posSlack[iRow] = -1;
     
    24742476                    int jRow = row[columnStart[iColumn]];
    24752477                    if (!columnLower[iColumn]) {
    2476                          if (element[columnStart[iColumn]] > 0.0 && posSlack[jRow] < 0)
    2477                               posSlack[jRow] = iColumn;
    2478                          else if (element[columnStart[iColumn]] < 0.0 && negSlack[jRow] < 0)
    2479                               negSlack[jRow] = iColumn;
    2480                     } else if (!columnUpper[iColumn]) {
     2478                      if (element[columnStart[iColumn]] > 0.0) {
     2479                        if(posSlack[jRow] < 0) {
     2480                          posSlack[jRow] = iColumn;
     2481                        } else {
     2482                          int jColumn=posSlack[jRow];
     2483                          while (nextPosSlack[jColumn]>=0)
     2484                            jColumn=nextPosSlack[jColumn];
     2485                          nextPosSlack[jColumn]=iColumn;
     2486                        }
     2487                      } else if (element[columnStart[iColumn]] < 0.0) {
     2488                        if(negSlack[jRow] < 0) {
     2489                          negSlack[jRow] = iColumn;
     2490                        } else {
     2491                          int jColumn=negSlack[jRow];
     2492                          while (nextNegSlack[jColumn]>=0)
     2493                            jColumn=nextNegSlack[jColumn];
     2494                          nextNegSlack[jColumn]=iColumn;
     2495                        }
     2496                      }
     2497                    } else if (!columnUpper[iColumn]&& false) {// out for testing
    24812498                         if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0)
    24822499                              posSlack[jRow] = iColumn;
     
    24872504          }
    24882505          // now see what that does to row solution
     2506          const double * objective = model2->objective();
    24892507          double * rowSolution = model2->primalRowSolution();
    24902508          CoinZeroN (rowSolution, numberRows);
     
    24982516                    int jColumn = posSlack[iRow];
    24992517                    if (jColumn >= 0) {
    2500                          if (columnSolution[jColumn])
     2518                      // sort if more than one
     2519                      int nPos=1;
     2520                      sort[0]=jColumn;
     2521                      weight[0]=objective[jColumn];
     2522                      while (nextPosSlack[jColumn]>=0) {
     2523                        jColumn=nextPosSlack[jColumn];
     2524                        sort[nPos]=jColumn;
     2525                        weight[nPos++]=objective[jColumn];
     2526                      }
     2527                      if (nPos>1) {
     2528                        CoinSort_2(weight,weight+nPos,sort);
     2529                        for (int i=0;i<nPos;i++) {
     2530                          double difference = lower[iRow] - rowSolution[iRow];
     2531                          jColumn=sort[i];
     2532                          double elementValue = element[columnStart[jColumn]];
     2533                          assert (elementValue > 0.0);
     2534                          double value=columnSolution[jColumn];
     2535                          double movement = CoinMin(difference / elementValue, columnUpper[jColumn]-value);
     2536                          columnSolution[jColumn] += movement;
     2537                          rowSolution[iRow] += movement * elementValue;
     2538                        }
     2539                        continue;
     2540                      }
     2541                        if (jColumn<0||columnSolution[jColumn])
    25012542                              continue;
    25022543                         double difference = lower[iRow] - rowSolution[iRow];
     
    25152556                    int jColumn = negSlack[iRow];
    25162557                    if (jColumn >= 0) {
    2517                          if (columnSolution[jColumn])
     2558                      // sort if more than one
     2559                      int nNeg=1;
     2560                      sort[0]=jColumn;
     2561                      weight[0]=objective[jColumn];
     2562                      while (nextNegSlack[jColumn]>=0) {
     2563                        jColumn=nextNegSlack[jColumn];
     2564                        sort[nNeg]=jColumn;
     2565                        weight[nNeg++]=objective[jColumn];
     2566                      }
     2567                      if (nNeg>1) {
     2568                        CoinSort_2(weight,weight+nNeg,sort);
     2569                        for (int i=0;i<nNeg;i++) {
     2570                          double difference = rowSolution[iRow]-upper[iRow];
     2571                          jColumn=sort[i];
     2572                          double elementValue = element[columnStart[jColumn]];
     2573                          assert (elementValue < 0.0);
     2574                          double value=columnSolution[jColumn];
     2575                          double movement = CoinMin(difference / -elementValue, columnUpper[jColumn]-value);
     2576                          columnSolution[jColumn] += movement;
     2577                          rowSolution[iRow] += movement * elementValue;
     2578                        }
     2579                        continue;
     2580                      }
     2581                        if (jColumn<0||columnSolution[jColumn])
    25182582                              continue;
    25192583                         double difference = upper[iRow] - rowSolution[iRow];
     
    27782842                         small.dual(0);
    27792843#endif
    2780                       if (small.problemStatus())
     2844                      if (small.problemStatus()) {
     2845                        int numberIterations=small.numberIterations();
    27812846                        small.dual(0);
     2847                        small.setNumberIterations(small.numberIterations()+numberIterations);
     2848                      }
    27822849#else
    27832850                         int numberColumns = small.numberColumns();
     
    28742941               if (iPass > 20)
    28752942                    sumArtificials = 0.0;
    2876                if ((small.objectiveValue()*optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 5 && sumArtificials < 1.0e-8) ||
     2943               if ((small.objectiveValue()*optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 5 && sumArtificials < 1.0e-8 && maxSprintPass<200) ||
    28772944                         (!small.numberIterations() && iPass) ||
    28782945                         iPass == maxSprintPass - 1 || small.status() == 3) {
     
    42414308                    } else {
    42424309                         // primal
    4243                          if (model_->infeasibilityCost() > 1.0e14)
    4244                               model_->setInfeasibilityCost(1.0e14);
     4310                         //if (model_->infeasibilityCost() > 1.0e14)
     4311                         //   model_->setInfeasibilityCost(1.0e14);
    42454312                         iSequence = out_[CLP_CYCLE-1];
    42464313                    }
  • trunk/Clp/src/ClpSolve.hpp

    r2078 r2259  
    212212     /// Set whole group
    213213     inline int presolveActions() const {
    214           return independentOptions_[1] & 0xffff;
     214          return independentOptions_[1] & 0xffffff;
    215215     }
    216216     inline void setPresolveActions(int action) {
    217           independentOptions_[1]  = (independentOptions_[1] & 0xffff0000) | (action & 0xffff);
     217          independentOptions_[1]  = (independentOptions_[1] & 0xff000000) | (action & 0xffffff);
    218218     }
    219219     /// Largest column for substitution (normally 3)
Note: See TracChangeset for help on using the changeset viewer.