Changeset 2235


Ignore:
Timestamp:
Jan 3, 2017 10:44:15 AM (2 years ago)
Author:
forrest
Message:

stuff for vector matrix

Location:
trunk/Clp/src
Files:
28 edited

Legend:

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

    r2227 r2235  
    17051705     parameters[numberParameters-1].append("Mumps_dummy");
    17061706#endif
     1707#ifdef PARDISO_BARRIER
     1708     parameters[numberParameters-1].append("Pardiso");
     1709#else
     1710     parameters[numberParameters-1].append("Pardiso_dummy");
     1711#endif
    17071712     parameters[numberParameters-1].setLonghelp
    17081713     (
     
    22412246     parameters[numberParameters++] =
    22422247          CbcOrClpParam("exper!iment", "Whether to use testing features",
    2243                         -1, 200, CBC_PARAM_INT_EXPERIMENT, 0);
     2248                        -1, 200000, CBC_PARAM_INT_EXPERIMENT, 0);
    22442249     parameters[numberParameters-1].setLonghelp
    22452250     (
  • trunk/Clp/src/ClpEventHandler.hpp

    r2157 r2235  
    6363          startOfIterationInDual,
    6464          updateDualsInDual,
     65          beforeDeleteRim,
    6566          endOfCreateRim,
    6667          slightlyInfeasible,
     
    6869          moreMiniPresolve,
    6970          modifyMatrixInMiniPostsolve,
     71          beforeChooseIncoming,
     72          afterChooseIncoming,
     73          beforeCreateNonLinear,
     74          afterCreateNonLinear,
    7075          startOfCrossover, // in Idiot
    7176          noTheta // At end (because no pivot)
  • trunk/Clp/src/ClpFactorization.cpp

    r2149 r2235  
    55
    66#include "CoinPragma.hpp"
     7#if ABOCA_LITE
     8// 1 is not owner of abcState_
     9#define ABCSTATE_LITE 1
     10#endif
    711#include "ClpFactorization.hpp"
    812#ifndef SLIM_CLP
     
    10081012          return networkBasis_->updateColumnTranspose(regionSparse, regionSparse2);
    10091013#endif
     1014     }
     1015#endif
     1016}
     1017/* Updates two columns (BTRAN) from regionSparse2 and 3
     1018   regionSparse starts as zero and is zero at end
     1019   Note - if regionSparse2 packed on input - will be packed on output - same for 3
     1020*/
     1021void
     1022ClpFactorization::updateTwoColumnsTranspose ( CoinIndexedVector * regionSparse,
     1023                                  CoinIndexedVector * regionSparse2,
     1024                                  CoinIndexedVector * regionSparse3) const
     1025{
     1026     if (!numberRows_)
     1027          return 0;
     1028#ifndef SLIM_CLP
     1029     if (!networkBasis_) {
     1030#endif
     1031#ifdef CLP_FACTORIZATION_INSTRUMENT
     1032          factorization_instrument(-1);
     1033#endif
     1034          collectStatistics_ = true;
     1035          CoinFactorization::updateTwoColumnsTranspose(regionSparse,
     1036                                                       regionSparse2,regionSparse3);
     1037          collectStatistics_ = false;
     1038#ifdef CLP_FACTORIZATION_INSTRUMENT
     1039          factorization_instrument(6);
     1040#endif
     1041          return returnCode;
     1042#ifndef SLIM_CLP
     1043     } else {
     1044       factorization__->updateColumnTranspose(regionSparse, regionSparse2);
     1045       factorization__->updateColumnTranspose(regionSparse, regionSparse3);
    10101046     }
    10111047#endif
     
    29683004          return 0;
    29693005     coinFactorizationA_->setCollectStatistics(false);
     3006     // above doesn't work any more - do by hand
     3007     double save[15];
     3008     memcpy(save,&coinFactorizationA_->ftranCountInput_,sizeof(save));
    29703009     int returnCode = coinFactorizationA_->updateColumn(regionSparse,
    29713010                      regionSparse2,
    29723011                      noPermute);
     3012     memcpy(&coinFactorizationA_->ftranCountInput_,save,sizeof(save));
    29733013     return returnCode;
    29743014}
     
    30483088#endif
    30493089}
     3090/* Updates two columns (BTRAN) from regionSparse2 and 3
     3091   regionSparse starts as zero and is zero at end
     3092   Note - if regionSparse2 packed on input - will be packed on output - same for 3
     3093*/
     3094void
     3095ClpFactorization::updateTwoColumnsTranspose ( CoinIndexedVector * regionSparse,
     3096                                  CoinIndexedVector * regionSparse2,
     3097                                  CoinIndexedVector * regionSparse3) const
     3098{
     3099  if (!numberRows())
     3100          return;
     3101#ifndef SLIM_CLP
     3102     if (!networkBasis_) {
     3103#endif
     3104#ifdef CLP_FACTORIZATION_INSTRUMENT
     3105          factorization_instrument(-1);
     3106#endif
     3107          if (coinFactorizationA_) {
     3108               coinFactorizationA_->setCollectStatistics(doStatistics_);
     3109#if ABOCA_LITE_FACTORIZATION
     3110               coinFactorizationA_->updateTwoColumnsTranspose(regionSparse,
     3111                                                              regionSparse2,regionSparse3,abcState());
     3112#else
     3113               coinFactorizationA_->updateTwoColumnsTranspose(regionSparse,
     3114              regionSparse2,regionSparse3,0);
     3115#endif
     3116               coinFactorizationA_->setCollectStatistics(false);
     3117          } else {
     3118               coinFactorizationB_->updateColumnTranspose(regionSparse,
     3119                            regionSparse2);
     3120               coinFactorizationB_->updateColumnTranspose(regionSparse,
     3121                            regionSparse3);
     3122          }
     3123#ifdef CLP_FACTORIZATION_INSTRUMENT
     3124          factorization_instrument(6);
     3125#endif
     3126#ifndef SLIM_CLP
     3127     } else {
     3128       updateColumnTranspose(regionSparse, regionSparse2);
     3129       updateColumnTranspose(regionSparse, regionSparse3);
     3130     }
     3131#endif
     3132}
    30503133/* makes a row copy of L for speed and to allow very sparse problems */
    30513134void
     
    31603243#endif
    31613244}
     3245#if ABOCA_LITE_FACTORIZATION
     3246// Does btranU part of replaceColumn (skipping entries)
     3247void
     3248ClpFactorization::replaceColumn1 ( CoinIndexedVector * regionSparse,
     3249                                    int pivotRow)
     3250{
     3251  if (coinFactorizationA_)
     3252    coinFactorizationA_->replaceColumn1(regionSparse,pivotRow);
     3253}
     3254// Does replaceColumn - having already done btranU
     3255int
     3256ClpFactorization::replaceColumn2 ( CoinIndexedVector * regionSparse,
     3257                                 int pivotRow,
     3258                                    double pivotCheck)
     3259{
     3260  if (coinFactorizationA_)
     3261    return coinFactorizationA_->replaceColumn2(regionSparse,pivotRow,pivotCheck);
     3262  else
     3263    return 12345678;
     3264}
     3265#endif
    31623266// Set tolerances to safer of existing and given
    31633267void
  • trunk/Clp/src/ClpFactorization.hpp

    r2149 r2235  
    9696                         bool checkBeforeModifying = false,
    9797                         double acceptablePivot = 1.0e-8);
     98#if ABOCA_LITE_FACTORIZATION
     99  /// Does btranU part of replaceColumn (skipping entries)
     100  void replaceColumn1(CoinIndexedVector * regionSparse, int pivotRow);
     101  /// Does replaceColumn - having already done btranU
     102  int replaceColumn2 ( CoinIndexedVector * regionSparse,
     103                      int pivotRow,
     104                       double pivotCheck);
     105#endif
    98106     //@}
    99107
     
    129137     int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    130138                                 CoinIndexedVector * regionSparse2) const;
     139     /** Updates two columns (BTRAN) from regionSparse2 and 3
     140         regionSparse starts as zero and is zero at end
     141         Note - if regionSparse2 packed on input - will be packed on output - same for 3
     142     */
     143     void updateTwoColumnsTranspose ( CoinIndexedVector * regionSparse,
     144                                      CoinIndexedVector * regionSparse2,
     145                                      CoinIndexedVector * regionSparse3) const;
    131146     //@}
    132147#ifdef CLP_MULTIPLE_FACTORIZATIONS
     
    348363     inline int isDenseOrSmall() const {
    349364          return coinFactorizationB_ ? 1 : 0;
     365     }
     366     /// Return coinFactorizationA_
     367     inline CoinFactorization * coinFactorization() const {
     368          return coinFactorizationA_;
    350369     }
    351370#else
  • trunk/Clp/src/ClpMatrixBase.cpp

    r1665 r2235  
    564564}
    565565// Updates two arrays for steepest
    566 void
     566int
    567567ClpMatrixBase::transposeTimes2(const ClpSimplex * ,
    568568                               const CoinIndexedVector * , CoinIndexedVector *,
    569569                               const CoinIndexedVector * ,
    570570                               CoinIndexedVector * ,
     571                               double * , double * ,
    571572                               double , double ,
    572573                               // Array for exact devex to say what is in reference framework
     
    576577     std::cerr << "transposeTimes2 not supported - ClpMatrixBase" << std::endl;
    577578     abort();
     579     return 0;
    578580}
    579581/* Set the dimensions of the matrix. In effect, append new empty
  • trunk/Clp/src/ClpMatrixBase.hpp

    r2078 r2235  
    111111         default does not allow scaling
    112112         returns non-zero if no scaling done */
    113      virtual int scale(ClpModel * , const ClpSimplex * = NULL) const {
     113     virtual int scale(ClpModel * , ClpSimplex * = NULL) const {
    114114          return 1;
    115115     }
     
    322322          return false;
    323323     }
    324      /// Updates two arrays for steepest and does devex weights (need not be coded)
    325      virtual void transposeTimes2(const ClpSimplex * model,
    326                                   const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    327                                   const CoinIndexedVector * pi2,
    328                                   CoinIndexedVector * spare,
     324     /** Updates two arrays for steepest and does devex weights
     325         (need not be coded)
     326         Returns nonzero if updates reduced cost and infeas -
     327         new infeas in dj1 */
     328     virtual int transposeTimes2(const ClpSimplex * model,
     329                                 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
     330                                 const CoinIndexedVector * pi2,
     331                                 CoinIndexedVector * spare,
     332                                 double * infeas, double * reducedCost,
    329333                                  double referenceIn, double devex,
    330334                                  // Array for exact devex to say what is in reference framework
  • trunk/Clp/src/ClpMessage.cpp

    r1926 r2235  
    3434     {CLP_DUAL_CHECKB, 12, 2, "New dual bound of %g"},
    3535     {CLP_DUAL_ORIGINAL, 13, 3, "Going back to original objective"},
    36      {CLP_SIMPLEX_PERTURB, 14, 1, "Perturbing problem by %g %% of %g - largest nonzero change %g (%% %g) - largest zero change %g"},
     36     {CLP_SIMPLEX_PERTURB, 14, 1, "Perturbing problem by %g%%of %g - largest nonzero change %g ( %g%%) - largest zero change %g"},
    3737     {CLP_PRIMAL_ORIGINAL, 15, 2, "Going back to original tolerance"},
    3838     {CLP_PRIMAL_WEIGHT, 16, 2, "New infeasibility weight of %g"},
  • trunk/Clp/src/ClpNonLinearCost.cpp

    r2030 r2235  
    1313#include "CoinHelperFunctions.hpp"
    1414#include "ClpNonLinearCost.hpp"
     15#include "ClpEventHandler.hpp"
    1516//#############################################################################
    1617// Constructors / Destructor / Assignment
     
    279280}
    280281#endif
     282// Refresh one- assuming regions OK
     283void
     284ClpNonLinearCost::refresh(int iSequence)
     285{
     286     double infeasibilityCost = model_->infeasibilityCost();
     287     double primalTolerance = model_->currentPrimalTolerance();
     288     double * cost = model_->costRegion();
     289     double * upper = model_->upperRegion();
     290     double * lower = model_->lowerRegion();
     291     double * solution = model_->solutionRegion();
     292     cost2_[iSequence] = cost[iSequence];
     293     double value = solution[iSequence];
     294     double lowerValue = lower[iSequence];
     295     double upperValue = upper[iSequence];
     296     if (value - upperValue <= primalTolerance) {
     297       if (value - lowerValue >= -primalTolerance) {
     298         // feasible
     299         status_[iSequence] = static_cast<unsigned char>(CLP_FEASIBLE | (CLP_SAME << 4));
     300         bound_[iSequence] = 0.0;
     301       } else {
     302         // below
     303         cost[iSequence] -= infeasibilityCost;
     304         status_[iSequence] = static_cast<unsigned char>(CLP_BELOW_LOWER | (CLP_SAME << 4));
     305         bound_[iSequence] = upperValue;
     306         upper[iSequence] = lowerValue;
     307         lower[iSequence] = -COIN_DBL_MAX;
     308       }
     309     } else {
     310       // above
     311       cost[iSequence] += infeasibilityCost;
     312       status_[iSequence] = static_cast<unsigned char>(CLP_ABOVE_UPPER | (CLP_SAME << 4));
     313       bound_[iSequence] = lowerValue;
     314       lower[iSequence] = upperValue;
     315       upper[iSequence] = COIN_DBL_MAX;
     316     }
     317     
     318}
    281319// Refreshes costs always makes row costs zero
    282320void
     
    16041642          double costValue = cost2_[iSequence];
    16051643          int iWhere = originalStatus(iStatus);
     1644#undef CLP_USER_DRIVEN
     1645#ifdef CLP_USER_DRIVEN
     1646          clpUserStruct info;
     1647          info.type=3;
     1648          info.sequence=iSequence;
     1649          info.value=value;
     1650          info.change=0;
     1651          info.lower=lowerValue;
     1652          info.upper=upperValue;
     1653          info.cost=costValue;
     1654#endif             
    16061655          if (iWhere == CLP_BELOW_LOWER) {
    16071656               lowerValue = upperValue;
     
    16761725               break;
    16771726          }
     1727#ifdef CLP_USER_DRIVEN
     1728          model_->eventHandler()->eventWithInfo(ClpEventHandler::pivotRow,
     1729                                                           &info);
     1730          if (info.change) {
     1731          }
     1732#endif             
    16781733     }
    16791734     changeCost_ += value * difference;
     
    17261781     double difference = 0.0;
    17271782     int direction = 0;
     1783#ifdef CLP_USER_DRIVEN
     1784     double saveLower = model_->lowerRegion()[iSequence];
     1785     double saveUpper = model_->upperRegion()[iSequence];
     1786     double saveCost = model_->costRegion()[iSequence];
     1787#endif
    17281788     if (CLP_METHOD1) {
    17291789          // get where in bound sequence
     
    18821942     }
    18831943     changeCost_ += value * difference;
     1944#ifdef CLP_USER_DRIVEN
     1945     //if (iSequence>=model_->numberColumns)
     1946     if (difference||saveCost!=model_->costRegion()[iSequence]) {
     1947       printf("Out %d (%d) %g <= %g <= %g cost %g x=>x %g < %g cost %g\n",
     1948              iSequence,iSequence-model_->numberColumns(),
     1949              saveLower,value,saveUpper,saveCost,
     1950              model_->lowerRegion()[iSequence],
     1951              model_->upperRegion()[iSequence],
     1952              model_->costRegion()[iSequence]);
     1953     }
     1954#endif
    18841955     return direction;
    18851956}
  • trunk/Clp/src/ClpNonLinearCost.hpp

    r1769 r2235  
    152152     /// Refresh - assuming regions OK
    153153     void refresh();
     154     /// Refresh one- assuming regions OK
     155     void refresh(int iSequence);
    154156     /** Sets bounds and cost for one variable
    155157         Returns change in cost
     
    282284          return cost_[whichRange_[sequence] + offset_[sequence]];
    283285     }
     286     /// Returns full status
     287     inline int fullStatus(int sequence) const {
     288       return status_[sequence];
     289     }
     290     /// Returns if changed from beginning of iteration
     291     inline bool changed(int sequence) const {
     292       return (status_[sequence]&64)==0;
     293     }
     294                   
    284295     //@}
    285296
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r2078 r2235  
    1212#include "CoinHelperFunctions.hpp"
    1313//#define THREAD
    14 
     14//#define FAKE_CILK
     15#if ABOCA_LITE
     16// 1 is not owner of abcState_
     17#define ABCSTATE_LITE 1
     18#endif
    1519#include "ClpSimplex.hpp"
    1620#include "ClpSimplexDual.hpp"
     
    4852// Constructors / Destructor / Assignment
    4953//#############################################################################
    50 
     54static void debug3(int iColumn,double &thisWeight,double pivotSquared,double devex,
     55                   double pivot,double modification,double oldWeight)
     56{
     57  //double newWeight=oldWeight+pivotSquared * devex + pivot * modification;
     58  //if (/*iColumn==34845 && */newWeight!=thisWeight)
     59  //printf("YY %d %.30g %.30g %.30g %.30g %.30g %.30g\n",
     60  //       iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight);
     61  thisWeight=oldWeight;
     62  thisWeight += pivotSquared * devex + pivot * modification;
     63}
    5164//-------------------------------------------------------------------
    5265// Default Constructor
     
    316329     }
    317330}
     331#if ABOCA_LITE
     332static void
     333transposeTimesBit(clpTempInfo & info)
     334{
     335  const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
     336  const int * COIN_RESTRICT row = info.row;
     337  const double * COIN_RESTRICT elementByColumn = info.element;
     338  double * COIN_RESTRICT y = info.spare;
     339  const double * COIN_RESTRICT x = info.work;
     340  int first = info.startColumn;
     341  int last = first+info.numberToDo;
     342  CoinBigIndex start = columnStart[first];
     343  for (int iColumn = first; iColumn < last; iColumn++) {
     344    CoinBigIndex j;
     345    CoinBigIndex next = columnStart[iColumn+1];
     346    double value = y[iColumn];
     347    for (j = start; j < next; j++) {
     348      int jRow = row[j];
     349      value -= x[jRow] * elementByColumn[j];
     350    }
     351    start = next;
     352    y[iColumn] = value;
     353  }
     354}
     355#endif
    318356void
    319357ClpPackedMatrix::transposeTimes(double scalar,
     
    327365     if (!(flags_ & 2)) {
    328366          if (scalar == -1.0) {
     367#if ABOCA_LITE
     368            int numberThreads=abcState();
     369            if (numberThreads) {
     370            clpTempInfo info[ABOCA_LITE];
     371            int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
     372            int n=0;
     373            for (int i=0;i<numberThreads;i++) {
     374              info[i].spare=y;
     375              info[i].work=const_cast<double *>(x);
     376              info[i].startColumn=n;
     377              info[i].element=elementByColumn;
     378              info[i].start=columnStart;
     379              info[i].row=row;
     380              info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
     381              n += chunk;
     382            }
     383            for (int i=0;i<numberThreads;i++) {
     384              cilk_spawn transposeTimesBit(info[i]);
     385            }
     386            cilk_sync;
     387            } else {
     388#endif
    329389               CoinBigIndex start = columnStart[0];
    330390               for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     
    339399                    y[iColumn] = value;
    340400               }
     401#if ABOCA_LITE
     402            }
     403#endif
    341404          } else {
    342405               CoinBigIndex start = columnStart[0];
     
    519582     }
    520583}
     584#if ABOCA_LITE
     585static void
     586transposeTimesSubsetBit(clpTempInfo & info)
     587{
     588  const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
     589  const int * COIN_RESTRICT row = info.row;
     590  const double * COIN_RESTRICT elementByColumn = info.element;
     591  double * COIN_RESTRICT y = info.spare;
     592  const double * COIN_RESTRICT x = info.work;
     593  const int * COIN_RESTRICT which = info.which;
     594  int first = info.startColumn;
     595  int last = first+info.numberToDo;
     596  for (int i = first; i < last; i++) {
     597    int iColumn = which[i];
     598    CoinBigIndex j;
     599    CoinBigIndex start = columnStart[iColumn];
     600    CoinBigIndex next = columnStart[iColumn+1];
     601    double value = 0.0;
     602    for (j = start; j < next; j++) {
     603      int jRow = row[j];
     604      value += x[jRow] * elementByColumn[j];
     605    }
     606    start = next;
     607    y[iColumn] -= value;
     608  }
     609}
     610#endif
    521611void
    522612ClpPackedMatrix::transposeTimesSubset( int number,
     
    546636               }
    547637          } else {
     638#if ABOCA_LITE
     639               int numberThreads=abcState();
     640               if (numberThreads) {
     641               clpTempInfo info[ABOCA_LITE];
     642               int chunk = (number+numberThreads-1)/numberThreads;
     643               int n=0;
     644               for (int i=0;i<numberThreads;i++) {
     645                 info[i].spare=y;
     646                 info[i].work=const_cast<double *>(x);
     647                 info[i].which=const_cast<int *>(which);
     648                 info[i].startColumn=n;
     649                 info[i].element=elementByColumn;
     650                 info[i].start=columnStart;
     651                 info[i].row=row;
     652                 info[i].numberToDo=CoinMin(chunk,number-n);
     653                 n += chunk;
     654               }
     655               for (int i=0;i<numberThreads;i++) {
     656                 cilk_spawn transposeTimesSubsetBit(info[i]);
     657               }
     658               cilk_sync;
     659               } else {
     660#endif
    548661               for (int jColumn = 0; jColumn < number; jColumn++) {
    549662                    int iColumn = which[jColumn];
     
    558671                    y[iColumn] -= value;
    559672               }
     673#if ABOCA_LITE
     674            }
     675#endif
    560676          }
    561677     } else {
     
    9631079                              model->spareDoubleArray_[1] = bestPossible;
    9641080                              spareArray->setNumElements(numberRemaining);
     1081#if 0
     1082                              columnArray->setNumElements(numberNonZero);
     1083                              spareArray->setPackedMode(true);
     1084                              columnArray->setPackedMode(true);
     1085                              spareArray->checkClean();
     1086                              columnArray->checkClean();
     1087                              printf("col %d spare %d upper %g best %g\n",
     1088                                     columnArray->getNumElements(),
     1089                                     spareArray->getNumElements(),upperTheta,bestPossible);
     1090#endif
    9651091                              // signal partially done
    9661092                              model->spareIntArray_[0] = -2;
     
    14561582     return numberNonZero;
    14571583}
     1584#if ABOCA_LITE
     1585static void
     1586transposeTimesUnscaledBit(clpTempInfo & info)
     1587{
     1588  const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
     1589  const int * COIN_RESTRICT row = info.row;
     1590  const double * COIN_RESTRICT elementByColumn = info.element;
     1591  double * COIN_RESTRICT array = info.infeas;
     1592  int * COIN_RESTRICT index = info.which;
     1593  const unsigned char * COIN_RESTRICT status = info.status;
     1594  double zeroTolerance=info.tolerance;
     1595  const double * COIN_RESTRICT pi = info.work;
     1596  int first = info.startColumn;
     1597  int last = first+info.numberToDo;
     1598  double value = 0.0;
     1599  int jColumn = -1;
     1600  int numberNonZero=0;
     1601  for (int iColumn = first; iColumn < last; iColumn++) {
     1602    bool wanted = ((status[iColumn] & 3) != 1);
     1603    if (fabs(value) > zeroTolerance) {
     1604      array[numberNonZero] = value;
     1605      index[numberNonZero++] = jColumn;
     1606    }
     1607    value = 0.0;
     1608    if (wanted) {
     1609      CoinBigIndex start = columnStart[iColumn];
     1610      CoinBigIndex end = columnStart[iColumn+1];
     1611      jColumn = iColumn;
     1612      int n = end - start;
     1613      bool odd = (n & 1) != 0;
     1614      n = n >> 1;
     1615      const int * COIN_RESTRICT rowThis = row + start;
     1616      const double * COIN_RESTRICT elementThis = elementByColumn + start;
     1617      for (; n; n--) {
     1618        int iRow0 = *rowThis;
     1619        int iRow1 = *(rowThis + 1);
     1620        rowThis += 2;
     1621        value += pi[iRow0] * (*elementThis);
     1622        value += pi[iRow1] * (*(elementThis + 1));
     1623        elementThis += 2;
     1624      }
     1625      if (odd) {
     1626        int iRow = *rowThis;
     1627        value += pi[iRow] * (*elementThis);
     1628      }
     1629    }
     1630  }
     1631  if (fabs(value) > zeroTolerance) {
     1632    array[numberNonZero] = value;
     1633    index[numberNonZero++] = jColumn;
     1634  }
     1635  info.numberAdded=numberNonZero;
     1636}
     1637#endif
    14581638// Meat of transposeTimes by column when not scaled
    14591639int
     
    14691649     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    14701650     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1651#if ABOCA_LITE
     1652     int numberThreads=abcState();
     1653     if (numberThreads) {
     1654     clpTempInfo info[ABOCA_LITE];
     1655     int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
     1656     int n=0;
     1657     for (int i=0;i<numberThreads;i++) {
     1658       info[i].which = index+n;
     1659       info[i].infeas = array+n;
     1660       info[i].work=const_cast<double *>(pi);
     1661       info[i].status=const_cast<unsigned char *>(status);
     1662       info[i].startColumn=n;
     1663       info[i].element=elementByColumn;
     1664       info[i].start=columnStart;
     1665       info[i].row=row;
     1666       info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
     1667       info[i].tolerance=zeroTolerance;
     1668       n += chunk;
     1669     }
     1670     for (int i=0;i<numberThreads;i++) {
     1671       cilk_spawn transposeTimesUnscaledBit(info[i]);
     1672     }
     1673     cilk_sync;
     1674     for (int i=0;i<numberThreads;i++)
     1675       numberNonZero += info[i].numberAdded;
     1676     moveAndZero(info,2,NULL);
     1677       } else {
     1678#endif
    14711679     double value = 0.0;
    14721680     int jColumn = -1;
     
    15051713          index[numberNonZero++] = jColumn;
    15061714     }
     1715#if ABOCA_LITE
     1716            }
     1717#endif
    15071718     return numberNonZero;
    15081719}
     1720#if ABOCA_LITE
     1721static void
     1722transposeTimesUnscaledBit2(clpTempInfo & info)
     1723{
     1724  const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
     1725  const int * COIN_RESTRICT row = info.row;
     1726  const double * COIN_RESTRICT elementByColumn = info.element;
     1727  double * COIN_RESTRICT array = info.infeas;
     1728  const double * COIN_RESTRICT reducedCost = info.reducedCost;
     1729  int * COIN_RESTRICT index = info.which;
     1730  double * COIN_RESTRICT spareArray = info.spare;
     1731  int * COIN_RESTRICT spareIndex = info.index;
     1732  const unsigned char * COIN_RESTRICT status = info.status;
     1733  double zeroTolerance=info.tolerance;
     1734  double dualTolerance=info.dualTolerance;
     1735  double acceptablePivot = info.acceptablePivot;
     1736  const double * COIN_RESTRICT pi = info.work;
     1737  double tentativeTheta = 1.0e15;
     1738  int numberRemaining = 0;
     1739  double upperTheta = info.upperTheta;
     1740  double bestPossible = info.bestPossible;
     1741  int first = info.startColumn;
     1742  int last = first+info.numberToDo;
     1743  int numberNonZero=0;
     1744  double multiplier[] = { -1.0, 1.0};
     1745  double dualT = - dualTolerance;
     1746  for (int iColumn = first; iColumn < last; iColumn++) {
     1747    int wanted = (status[iColumn] & 3) - 1;
     1748    if (wanted) {
     1749      double value = 0.0;
     1750      CoinBigIndex start = columnStart[iColumn];
     1751      CoinBigIndex end = columnStart[iColumn+1];
     1752      int n = end - start;
     1753      bool odd = (n & 1) != 0;
     1754      n = n >> 1;
     1755      const int * COIN_RESTRICT rowThis = row + start;
     1756      const double * COIN_RESTRICT elementThis = elementByColumn + start;
     1757      for (; n; n--) {
     1758        int iRow0 = *rowThis;
     1759        int iRow1 = *(rowThis + 1);
     1760        rowThis += 2;
     1761        value += pi[iRow0] * (*elementThis);
     1762        value += pi[iRow1] * (*(elementThis + 1));
     1763        elementThis += 2;
     1764      }
     1765      if (odd) {
     1766        int iRow = *rowThis;
     1767        value += pi[iRow] * (*elementThis);
     1768      }
     1769      if (fabs(value) > zeroTolerance) {
     1770        double mult = multiplier[wanted-1];
     1771        double alpha = value * mult;
     1772        array[numberNonZero] = value;
     1773        index[numberNonZero++] = iColumn;
     1774        if (alpha > 0.0) {
     1775          double oldValue = reducedCost[iColumn] * mult;
     1776          double value = oldValue - tentativeTheta * alpha;
     1777          if (value < dualT) {
     1778            bestPossible = CoinMax(bestPossible, alpha);
     1779            value = oldValue - upperTheta * alpha;
     1780            if (value < dualT && alpha >= acceptablePivot) {
     1781              upperTheta = (oldValue - dualT) / alpha;
     1782              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     1783            }
     1784            // add to list
     1785            spareArray[numberRemaining] = alpha * mult;
     1786            spareIndex[numberRemaining++] = iColumn;
     1787          }
     1788        }
     1789      }
     1790    }
     1791  }
     1792  info.numberAdded=numberNonZero;
     1793  info.numberRemaining=numberRemaining;
     1794  info.bestPossible=bestPossible;
     1795  info.upperTheta=upperTheta;
     1796}
     1797#endif
    15091798/* Meat of transposeTimes by column when not scaled and skipping
    15101799   and doing part of dualColumn */
     
    15241813          const double zeroTolerance) const
    15251814{
    1526      double tentativeTheta = 1.0e15;
    15271815     int numberRemaining = numberRemainingP;
    15281816     double upperTheta = upperThetaP;
     
    15331821     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
    15341822     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1823#if ABOCA_LITE
     1824     int numberThreads=abcState();
     1825     if (numberThreads) {
     1826     clpTempInfo info[ABOCA_LITE];
     1827     int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
     1828     int n=0;
     1829     for (int i=0;i<numberThreads;i++) {
     1830       info[i].upperTheta = upperThetaP;
     1831       info[i].bestPossible = bestPossibleP;
     1832       info[i].acceptablePivot = acceptablePivot;
     1833       info[i].reducedCost = const_cast<double *>(reducedCost);
     1834       info[i].index = spareIndex+n+numberRemaining;
     1835       info[i].spare = spareArray+n+numberRemaining;
     1836       info[i].which = index+n;
     1837       info[i].infeas = array+n;
     1838       info[i].work=const_cast<double *>(pi);
     1839       info[i].status=const_cast<unsigned char *>(status);
     1840       info[i].startColumn=n;
     1841       info[i].element=elementByColumn;
     1842       info[i].start=columnStart;
     1843       info[i].row=row;
     1844       info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
     1845       info[i].tolerance=zeroTolerance;
     1846       info[i].dualTolerance=dualTolerance;
     1847       n += chunk;
     1848     }
     1849     for (int i=0;i<numberThreads;i++) {
     1850       cilk_spawn transposeTimesUnscaledBit2(info[i]);
     1851     }
     1852     cilk_sync;
     1853     for (int i=0;i<numberThreads;i++) {
     1854       numberNonZero += info[i].numberAdded;
     1855       numberRemaining += info[i].numberRemaining;
     1856       bestPossible = CoinMax(bestPossible,static_cast<double>(info[i].bestPossible));
     1857       upperTheta = CoinMin(upperTheta,static_cast<double>(info[i].upperTheta));
     1858     }
     1859     moveAndZero(info,1,NULL);
     1860     moveAndZero(info,2,NULL);
     1861     } else {
     1862#endif
     1863     double tentativeTheta = 1.0e15;
    15351864     double multiplier[] = { -1.0, 1.0};
    15361865     double dualT = - dualTolerance;
     
    16591988          }
    16601989     }
     1990#if ABOCA_LITE
     1991            }
     1992#endif
    16611993     numberRemainingP = numberRemaining;
    16621994     upperThetaP = upperTheta;
     
    17042036     return numberNonZero;
    17052037}
     2038#if ABOCA_LITE
     2039static void
     2040packDownBit(clpTempInfo & info)
     2041{
     2042  double * COIN_RESTRICT output = info.infeas;
     2043  int * COIN_RESTRICT index = info.which;
     2044  double zeroTolerance=info.tolerance;
     2045  int first = info.startColumn;
     2046  int last = info.numberToDo;
     2047  int numberNonZero=0;
     2048  for (int i = 0; i < last; i++) {
     2049    double value = output[i];
     2050    if (value) {
     2051      output[i] = 0.0;
     2052      if (fabs(value) > zeroTolerance) {
     2053        output[numberNonZero] = value;
     2054        index[numberNonZero++] = i+first;
     2055      }
     2056    }
     2057  }
     2058  info.numberAdded=numberNonZero;
     2059}
     2060#endif
    17062061// Meat of transposeTimes by row n > K if packed - returns number nonzero
    17072062int
     
    17412096     // get rid of tiny values and count
    17422097     int numberNonZero = 0;
     2098#if ABOCA_LITE
     2099     int numberThreads=abcState();
     2100     if (numberThreads) {
     2101     clpTempInfo info[ABOCA_LITE];
     2102     int chunk = (numberColumns+numberThreads-1)/numberThreads;
     2103     int n=0;
     2104     for (int i=0;i<numberThreads;i++) {
     2105       info[i].which = index+n;
     2106       info[i].infeas = output+n;
     2107       info[i].startColumn=n;
     2108       info[i].numberToDo=CoinMin(chunk,numberColumns-n);
     2109       info[i].tolerance=tolerance;
     2110       n += chunk;
     2111     }
     2112     for (int i=0;i<numberThreads;i++) {
     2113       cilk_spawn packDownBit(info[i]);
     2114     }
     2115     cilk_sync;
     2116     for (int i=0;i<numberThreads;i++) {
     2117       numberNonZero += info[i].numberAdded;
     2118     }
     2119     moveAndZero(info,2,NULL);
     2120     } else {
     2121#endif
    17432122     for (int i = 0; i < numberColumns; i++) {
    17442123          double value = output[i];
     
    17512130          }
    17522131     }
     2132#if ABOCA_LITE
     2133            }
     2134#endif
    17532135#ifndef NDEBUG
    17542136     for (int i = numberNonZero; i < numberColumns; i++)
     
    20492431     if (!packed)
    20502432          factor *= 0.9;
     2433     // bias if columnCopy
     2434     if (columnCopy_)
     2435       factor *= 0.5;
    20512436     return ((numberInRowArray > factor * numberRows || !model->rowCopy()) && !(flags_ & 2));
    20522437}
     
    20552440#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
    20562441#endif
     2442#if ABOCA_LITE
     2443static void
     2444transposeTimes2UnscaledBit(clpTempInfo & info)
     2445{
     2446  const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
     2447  const int * COIN_RESTRICT row = info.row;
     2448  const double * COIN_RESTRICT elementByColumn = info.element;
     2449  double * COIN_RESTRICT array = info.infeas;
     2450  double * COIN_RESTRICT weights = const_cast<double *>(info.lower);
     2451  double * COIN_RESTRICT reducedCost = info.reducedCost;
     2452  double * COIN_RESTRICT infeas = const_cast<double *>(info.upper);
     2453  int * COIN_RESTRICT index = info.which;
     2454  unsigned int * COIN_RESTRICT reference = reinterpret_cast<unsigned int *>(info.index);
     2455  const unsigned char * COIN_RESTRICT status = info.status;
     2456  double zeroTolerance=info.tolerance;
     2457  double dualTolerance=info.dualTolerance;
     2458  const double * COIN_RESTRICT pi = info.work;
     2459  const double * COIN_RESTRICT piWeight = info.cost;
     2460  double scaleFactor = info.primalRatio;
     2461  double devex = info.changeObj;
     2462  double referenceIn = info.theta;
     2463  int first = info.startColumn;
     2464  int last = first+info.numberToDo;
     2465  int killDjs=info.numberInfeasibilities;
     2466  int numberNonZero=0;
     2467  double mmult[2]={1.0,-1.0};
     2468  for (int iColumn = first; iColumn < last; iColumn++) {
     2469    int iStat = status[iColumn]&7;
     2470    if (iStat!=1&&iStat!=5) {
     2471      double value = 0.0;
     2472      CoinBigIndex start = columnStart[iColumn];
     2473      CoinBigIndex end = columnStart[iColumn+1];
     2474      int n = end - start;
     2475      const int * COIN_RESTRICT rowThis = row + start;
     2476      const double * COIN_RESTRICT elementThis = elementByColumn + start;
     2477      for (int i=0; i < n; i ++) {
     2478        int iRow = rowThis[i];
     2479        value -= pi[iRow] * elementThis[i];
     2480      }
     2481      if (fabs(value) > zeroTolerance) {
     2482        // and do other array
     2483        double modification = 0.0;
     2484        for (int i=0; i < n; i ++) {
     2485          int iRow = rowThis[i];
     2486          modification += piWeight[iRow] * elementThis[i];
     2487        }
     2488        double thisWeight = weights[iColumn];
     2489        double oldWeight=thisWeight;
     2490        double pivot = value * scaleFactor;
     2491        double pivotSquared = pivot * pivot;
     2492        thisWeight += pivotSquared * devex + pivot * modification;
     2493        debug3(iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight);
     2494        if (thisWeight < DEVEX_TRY_NORM) {
     2495          if (referenceIn < 0.0) {
     2496            // steepest
     2497            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2498          } else {
     2499            // exact
     2500            thisWeight = referenceIn * pivotSquared;
     2501            if (reference(iColumn))
     2502              thisWeight += 1.0;
     2503            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2504          }
     2505        }
     2506        weights[iColumn] = thisWeight;
     2507        if (!killDjs) {
     2508          value = reducedCost[iColumn]-value;
     2509          reducedCost[iColumn] = value;
     2510          if (iStat>1&&iStat<4) {
     2511            value *= mmult[iStat-2];
     2512            if (value<=dualTolerance)
     2513              value=0.0;
     2514          } else {
     2515            // free or super basic
     2516            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2517              // we are going to bias towards free (but only if reasonable)
     2518              value *= FREE_BIAS;
     2519            } else {
     2520              value=0.0;
     2521            }
     2522          }
     2523          if (value) {
     2524            value *= value;
     2525            // store square in list
     2526            if (infeas[iColumn]) {
     2527              infeas[iColumn] = value; // already there
     2528            } else {
     2529              array[numberNonZero]=value;
     2530              index[numberNonZero++]=iColumn;
     2531            }
     2532          } else {
     2533            array[numberNonZero]=0.0;
     2534            index[numberNonZero++]=iColumn;
     2535          }
     2536        }
     2537      }
     2538    }
     2539  }
     2540  info.numberAdded=numberNonZero;
     2541}
     2542static void
     2543transposeTimes2ScaledBit(clpTempInfo & info)
     2544{
     2545  const CoinBigIndex * COIN_RESTRICT columnStart=info.start;
     2546  const int * COIN_RESTRICT row = info.row;
     2547  const double * COIN_RESTRICT elementByColumn = info.element;
     2548  double * COIN_RESTRICT array = info.infeas;
     2549  double * COIN_RESTRICT weights = const_cast<double *>(info.lower);
     2550  int * COIN_RESTRICT index = info.which;
     2551  unsigned int * COIN_RESTRICT reference = reinterpret_cast<unsigned int *>(info.index);
     2552  const unsigned char * COIN_RESTRICT status = info.status;
     2553  double zeroTolerance=info.tolerance;
     2554  double dualTolerance=info.dualTolerance;
     2555  double * COIN_RESTRICT reducedCost = info.reducedCost;
     2556  double * COIN_RESTRICT infeas = const_cast<double *>(info.upper);
     2557  const double * COIN_RESTRICT pi = info.work;
     2558  const double * COIN_RESTRICT piWeight = info.cost;
     2559  const double * COIN_RESTRICT columnScale = info.solution;
     2560  double scaleFactor = info.primalRatio;
     2561  double devex = info.changeObj;
     2562  double referenceIn = info.theta;
     2563  int first = info.startColumn;
     2564  int last = first+info.numberToDo;
     2565  int killDjs=info.numberInfeasibilities;
     2566  int numberNonZero=0;
     2567  double mmult[2]={1.0,-1.0};
     2568  for (int iColumn = first; iColumn < last; iColumn++) {
     2569    int iStat = status[iColumn]&7;
     2570    if (iStat!=1&&iStat!=5) {
     2571      double value = 0.0;
     2572      double scale = columnScale[iColumn];
     2573      CoinBigIndex start = columnStart[iColumn];
     2574      CoinBigIndex end = columnStart[iColumn+1];
     2575      int n = end - start;
     2576      bool odd = (n & 1) != 0;
     2577      n &= ~1;
     2578      const int * COIN_RESTRICT rowThis = row + start;
     2579      const double * COIN_RESTRICT elementThis = elementByColumn + start;
     2580      for (int i=0; i < n; i += 2) {
     2581        int iRow0 = rowThis[i];
     2582        int iRow1 = rowThis[i+1];
     2583        value -= pi[iRow0] * elementThis[i];
     2584        value -= pi[iRow1] * elementThis[i+1];
     2585      }
     2586      if (odd) {
     2587        int iRow = rowThis[n];
     2588        value -= pi[iRow] * elementThis[n];
     2589      }
     2590      value *= scale;
     2591      if (fabs(value) > zeroTolerance) {
     2592        // and do other array
     2593        double modification = 0.0;
     2594        for (int i=0; i < n; i += 2) {
     2595          int iRow0 = rowThis[i];
     2596          int iRow1 = rowThis[i+1];
     2597          modification += piWeight[iRow0] * elementThis[i];
     2598          modification += piWeight[iRow1] * elementThis[i+1];
     2599        }
     2600        if (odd) {
     2601          int iRow = rowThis[n];
     2602          modification += piWeight[iRow] * elementThis[n];
     2603        }
     2604        modification *= scale;
     2605        double thisWeight = weights[iColumn];
     2606        double pivot = value * scaleFactor;
     2607        double pivotSquared = pivot * pivot;
     2608        thisWeight += pivotSquared * devex + pivot * modification;
     2609        if (thisWeight < DEVEX_TRY_NORM) {
     2610          if (referenceIn < 0.0) {
     2611            // steepest
     2612            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2613          } else {
     2614            // exact
     2615            thisWeight = referenceIn * pivotSquared;
     2616            if (reference(iColumn))
     2617              thisWeight += 1.0;
     2618            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2619          }
     2620        }
     2621        weights[iColumn] = thisWeight;
     2622        if (!killDjs) {
     2623          value = reducedCost[iColumn]-value;
     2624          reducedCost[iColumn] = value;
     2625          if (iStat>1&&iStat<4) {
     2626            value *= mmult[iStat-2];
     2627            if (value<=dualTolerance)
     2628              value=0.0;
     2629          } else {
     2630            // free or super basic
     2631            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2632              // we are going to bias towards free (but only if reasonable)
     2633              value *= FREE_BIAS;
     2634            } else {
     2635              value=0.0;
     2636            }
     2637          }
     2638          if (value) {
     2639            value *= value;
     2640            // store square in list
     2641            if (infeas[iColumn]) {
     2642              infeas[iColumn] = value; // already there
     2643            } else {
     2644              array[numberNonZero]=value;
     2645              index[numberNonZero++]=iColumn;
     2646            }
     2647          } else {
     2648            array[numberNonZero]=0.0;
     2649            index[numberNonZero++]=iColumn;
     2650          }
     2651        }
     2652      }
     2653    }
     2654  }
     2655  info.numberAdded=numberNonZero;
     2656}
     2657#endif
     2658/* Updates two arrays for steepest and does devex weights
     2659   Returns nonzero if updates reduced cost and infeas -
     2660   new infeas in dj1 */
    20572661// Updates two arrays for steepest
    2058 void
     2662int
    20592663ClpPackedMatrix::transposeTimes2(const ClpSimplex * model,
    20602664                                 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    20612665                                 const CoinIndexedVector * pi2,
    20622666                                 CoinIndexedVector * spare,
     2667                                 double * COIN_RESTRICT infeas,
     2668                                 double * COIN_RESTRICT reducedCost,
    20632669                                 double referenceIn, double devex,
    20642670                                 // Array for exact devex to say what is in reference framework
     
    20662672                                 double * COIN_RESTRICT weights, double scaleFactor)
    20672673{
     2674     int returnCode=0;
    20682675     // put row of tableau in dj1
    20692676     double * COIN_RESTRICT pi = pi1->denseVector();
     
    20732680     int numberInRowArray = pi1->getNumElements();
    20742681     double zeroTolerance = model->zeroTolerance();
     2682     double dualTolerance = model->currentDualTolerance();
     2683     // we can't really trust infeasibilities if there is dual error
     2684     // this coding has to mimic coding in checkDualSolution
     2685     double error = CoinMin(1.0e-2, model->largestDualError());
     2686     // allow tolerance at least slightly bigger than standard
     2687     dualTolerance = dualTolerance +  error;
    20752688     bool packed = pi1->packedMode();
    20762689     // do by column
     
    21092722                    pi[iRow] = piOld[i];
    21102723               }
    2111                if (!columnCopy_) {
     2724               if (!columnCopy_ || killDjs) {
     2725                 if (infeas)
     2726                   returnCode=1;
     2727#if ABOCA_LITE
     2728               int numberThreads=abcState();
     2729               if (numberThreads) {
     2730                 clpTempInfo info[ABOCA_LITE];
     2731                 int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
     2732                 int n=0;
     2733                 for (int i=0;i<numberThreads;i++) {
     2734                   info[i].theta = referenceIn;
     2735                   info[i].changeObj = devex;
     2736                   info[i].primalRatio = scaleFactor;
     2737                   info[i].lower = weights;
     2738                   info[i].reducedCost = reducedCost;
     2739                   info[i].upper = infeas;
     2740                   info[i].which = index+n;
     2741                   info[i].infeas = array+n;
     2742                   info[i].index = reinterpret_cast<int *>(reference);
     2743                   info[i].work=const_cast<double *>(pi);
     2744                   info[i].cost=const_cast<double *>(piWeight);
     2745                   info[i].status=const_cast<unsigned char *>(model->statusArray());
     2746                   info[i].startColumn=n;
     2747                   info[i].element=elementByColumn;
     2748                   info[i].start=columnStart;
     2749                   info[i].row=row;
     2750                   info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
     2751                   info[i].tolerance=zeroTolerance;
     2752                   info[i].dualTolerance=dualTolerance;
     2753                   info[i].numberInfeasibilities=killDjs ? 1 : 0;
     2754                   n += chunk;
     2755                 }
     2756                 for (int i=0;i<numberThreads;i++) {
     2757                   cilk_spawn transposeTimes2UnscaledBit(info[i]);
     2758                 }
     2759                 cilk_sync;
     2760                 for (int i=0;i<numberThreads;i++) {
     2761                   numberNonZero += info[i].numberAdded;
     2762                 }
     2763                 moveAndZero(info,2,NULL);
     2764               } else {
     2765#endif
    21122766                    CoinBigIndex j;
    21132767                    CoinBigIndex end = columnStart[0];
    21142768                    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    2115                          CoinBigIndex start = end;
    2116                          end = columnStart[iColumn+1];
    2117                          ClpSimplex::Status status = model->getStatus(iColumn);
    2118                          if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
    2119                          double value = 0.0;
    2120                          for (j = start; j < end; j++) {
    2121                               int iRow = row[j];
    2122                               value -= pi[iRow] * elementByColumn[j];
    2123                          }
    2124                          if (fabs(value) > zeroTolerance) {
    2125                               // and do other array
    2126                               double modification = 0.0;
    2127                               for (j = start; j < end; j++) {
    2128                                    int iRow = row[j];
    2129                                    modification += piWeight[iRow] * elementByColumn[j];
    2130                               }
    2131                               double thisWeight = weights[iColumn];
    2132                               double pivot = value * scaleFactor;
    2133                               double pivotSquared = pivot * pivot;
    2134                               thisWeight += pivotSquared * devex + pivot * modification;
    2135                               if (thisWeight < DEVEX_TRY_NORM) {
    2136                                    if (referenceIn < 0.0) {
    2137                                         // steepest
    2138                                         thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    2139                                    } else {
    2140                                         // exact
    2141                                         thisWeight = referenceIn * pivotSquared;
    2142                                         if (reference(iColumn))
    2143                                              thisWeight += 1.0;
    2144                                         thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    2145                                    }
    2146                               }
    2147                               weights[iColumn] = thisWeight;
    2148                               if (!killDjs) {
    2149                                    array[numberNonZero] = value;
    2150                                    index[numberNonZero++] = iColumn;
    2151                               }
    2152                          }
    2153                     }
     2769                      CoinBigIndex start = end;
     2770                      end = columnStart[iColumn+1];
     2771                      ClpSimplex::Status status = model->getStatus(iColumn);
     2772                      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
     2773                      double value = 0.0;
     2774                      for (j = start; j < end; j++) {
     2775                        int iRow = row[j];
     2776                        value -= pi[iRow] * elementByColumn[j];
     2777                      }
     2778                      if (fabs(value) > zeroTolerance) {
     2779                        // and do other array
     2780                        double modification = 0.0;
     2781                        for (j = start; j < end; j++) {
     2782                          int iRow = row[j];
     2783                          modification += piWeight[iRow] * elementByColumn[j];
     2784                        }
     2785                        double thisWeight = weights[iColumn];
     2786                        double oldWeight=thisWeight;
     2787                        double pivot = value * scaleFactor;
     2788                        double pivotSquared = pivot * pivot;
     2789                        thisWeight += pivotSquared * devex + pivot * modification;
     2790                        debug3(iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight);
     2791                        if (thisWeight < DEVEX_TRY_NORM) {
     2792                          if (referenceIn < 0.0) {
     2793                            // steepest
     2794                            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     2795                          } else {
     2796                            // exact
     2797                            thisWeight = referenceIn * pivotSquared;
     2798                            if (reference(iColumn))
     2799                              thisWeight += 1.0;
     2800                            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     2801                          }
     2802                        }
     2803                        weights[iColumn] = thisWeight;
     2804                        if (!killDjs) {
     2805                          value = reducedCost[iColumn]-value;
     2806                          reducedCost[iColumn] = value;
     2807                          // simplify status
     2808                          ClpSimplex::Status status = model->getStatus(iColumn);
     2809                         
     2810                          switch(status) {
     2811                           
     2812                          case ClpSimplex::basic:
     2813                          case ClpSimplex::isFixed:
     2814                            break;
     2815                          case ClpSimplex::isFree:
     2816                          case ClpSimplex::superBasic:
     2817                            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2818                              // we are going to bias towards free (but only if reasonable)
     2819                              value *= FREE_BIAS;
     2820                              value *= value;
     2821                              // store square in list
     2822                              if (infeas[iColumn]) {
     2823                                infeas[iColumn] = value; // already there
     2824                              } else {
     2825                                array[numberNonZero]=value;
     2826                                index[numberNonZero++]=iColumn;
     2827                              }
     2828                            } else {
     2829                              array[numberNonZero]=0.0;
     2830                              index[numberNonZero++]=iColumn;
     2831                            }
     2832                            break;
     2833                          case ClpSimplex::atUpperBound:
     2834                            if (value > dualTolerance) {
     2835                              value *= value;
     2836                              // store square in list
     2837                              if (infeas[iColumn]) {
     2838                                infeas[iColumn] = value; // already there
     2839                              } else {
     2840                                array[numberNonZero]=value;
     2841                                index[numberNonZero++]=iColumn;
     2842                              }
     2843                            } else {
     2844                              array[numberNonZero]=0.0;
     2845                              index[numberNonZero++]=iColumn;
     2846                            }
     2847                            break;
     2848                          case ClpSimplex::atLowerBound:
     2849                            if (value < -dualTolerance) {
     2850                              value *= value;
     2851                              // store square in list
     2852                              if (infeas[iColumn]) {
     2853                                infeas[iColumn] = value; // already there
     2854                              } else {
     2855                                array[numberNonZero]=value;
     2856                                index[numberNonZero++]=iColumn;
     2857                              }
     2858                            } else {
     2859                              array[numberNonZero]=0.0;
     2860                              index[numberNonZero++]=iColumn;
     2861                            }
     2862                          }
     2863                        }
     2864                      }
     2865                    }
     2866#if ABOCA_LITE
     2867            }
     2868#endif
    21542869               } else {
    21552870                    // use special column copy
    21562871                    // reset back
    2157                     if (killDjs)
    2158                          scaleFactor = 0.0;
    2159                     columnCopy_->transposeTimes2(model, pi, dj1, piWeight, referenceIn, devex,
     2872                    assert (!killDjs);
     2873                    //if (killDjs)
     2874                    //   scaleFactor = 0.0;
     2875                 if (infeas)
     2876                   returnCode=1;
     2877                    columnCopy_->transposeTimes2(model, pi, dj1, piWeight,
     2878                                                 infeas, reducedCost,
     2879                                                 referenceIn, devex,
    21602880                                                 reference, weights, scaleFactor);
    21612881                    numberNonZero = dj1->getNumElements();
     
    21752895                    piWeight[iRow] *= rowScale[iRow];
    21762896               }
    2177                if (!columnCopy_) {
    2178                     const double * COIN_RESTRICT columnScale = model->columnScale();
     2897               if (!columnCopy_ || killDjs) {
     2898                 if (infeas)
     2899                   returnCode=1;
     2900                 const double * COIN_RESTRICT columnScale = model->columnScale();
     2901#if ABOCA_LITE
     2902                 int numberThreads=abcState();
     2903                 if (numberThreads) {
     2904                   clpTempInfo info[ABOCA_LITE];
     2905                   int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads;
     2906                   int n=0;
     2907                   for (int i=0;i<numberThreads;i++) {
     2908                     info[i].theta = referenceIn;
     2909                     info[i].changeObj = devex;
     2910                     info[i].primalRatio = scaleFactor;
     2911                     info[i].lower = weights;
     2912                     info[i].reducedCost = reducedCost;
     2913                     info[i].upper = infeas;
     2914                     info[i].solution = const_cast<double *>(columnScale);
     2915                     info[i].which = index+n;
     2916                     info[i].infeas = array+n;
     2917                     info[i].index = reinterpret_cast<int *>(reference);
     2918                     info[i].work=const_cast<double *>(pi);
     2919                     info[i].cost=const_cast<double *>(piWeight);
     2920                     info[i].status=const_cast<unsigned char *>(model->statusArray());
     2921                     info[i].startColumn=n;
     2922                     info[i].element=elementByColumn;
     2923                     info[i].start=columnStart;
     2924                     info[i].row=row;
     2925                     info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n);
     2926                     info[i].tolerance=zeroTolerance;
     2927                     info[i].dualTolerance=dualTolerance;
     2928                     info[i].numberInfeasibilities=killDjs ? 1 : 0;
     2929                     n += chunk;
     2930                   }
     2931                   for (int i=0;i<numberThreads;i++) {
     2932                     cilk_spawn transposeTimes2ScaledBit(info[i]);
     2933                   }
     2934                   cilk_sync;
     2935                   for (int i=0;i<numberThreads;i++) {
     2936                     numberNonZero += info[i].numberAdded;
     2937                   }
     2938                   moveAndZero(info,2,NULL);
     2939                   if (infeas) {
     2940                     returnCode=1;
     2941                     dj1->setNumElements(numberNonZero);
     2942                   }
     2943                 } else {
     2944#endif
    21792945                    CoinBigIndex j;
    21802946                    CoinBigIndex end = columnStart[0];
     
    22162982                              weights[iColumn] = thisWeight;
    22172983                              if (!killDjs) {
    2218                                    array[numberNonZero] = value;
    2219                                    index[numberNonZero++] = iColumn;
    2220                               }
    2221                          }
    2222                     }
     2984                                value = reducedCost[iColumn]-value;
     2985                                reducedCost[iColumn] = value;
     2986                                // simplify status
     2987                                ClpSimplex::Status status = model->getStatus(iColumn);
     2988                               
     2989                                switch(status) {
     2990                                 
     2991                                case ClpSimplex::basic:
     2992                                case ClpSimplex::isFixed:
     2993                                  break;
     2994                                case ClpSimplex::isFree:
     2995                                case ClpSimplex::superBasic:
     2996                                  if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     2997                                    // we are going to bias towards free (but only if reasonable)
     2998                                    value *= FREE_BIAS;
     2999                                    value *= value;
     3000                                    // store square in list
     3001                                    if (infeas[iColumn]) {
     3002                                      infeas[iColumn] = value; // already there
     3003                                    } else {
     3004                                      array[numberNonZero]=value;
     3005                                      index[numberNonZero++]=iColumn;
     3006                                    }
     3007                                  } else {
     3008                                    array[numberNonZero]=0.0;
     3009                                    index[numberNonZero++]=iColumn;
     3010                                  }
     3011                                  break;
     3012                                case ClpSimplex::atUpperBound:
     3013                                  if (value > dualTolerance) {
     3014                                    value *= value;
     3015                                    // store square in list
     3016                                    if (infeas[iColumn]) {
     3017                                      infeas[iColumn] = value; // already there
     3018                                    } else {
     3019                                      array[numberNonZero]=value;
     3020                                      index[numberNonZero++]=iColumn;
     3021                                    }
     3022                                  } else {
     3023                                    array[numberNonZero]=0.0;
     3024                                    index[numberNonZero++]=iColumn;
     3025                                  }
     3026                                  break;
     3027                                case ClpSimplex::atLowerBound:
     3028                                  if (value < -dualTolerance) {
     3029                                    value *= value;
     3030                                    // store square in list
     3031                                    if (infeas[iColumn]) {
     3032                                      infeas[iColumn] = value; // already there
     3033                                    } else {
     3034                                      array[numberNonZero]=value;
     3035                                      index[numberNonZero++]=iColumn;
     3036                                    }
     3037                                  } else {
     3038                                    array[numberNonZero]=0.0;
     3039                                    index[numberNonZero++]=iColumn;
     3040                                  }
     3041                                }
     3042                              }
     3043                         }
     3044                    }
     3045#if ABOCA_LITE
     3046                 }
     3047#endif
     3048                 if (infeas&&false) {
     3049                   double tolerance = model->currentDualTolerance();
     3050                   // we can't really trust infeasibilities if there is dual error
     3051                   // this coding has to mimic coding in checkDualSolution
     3052                   double error = CoinMin(1.0e-2, model->largestDualError());
     3053                   // allow tolerance at least slightly bigger than standard
     3054                   tolerance = tolerance +  error;
     3055                   returnCode=1;
     3056                   int number=0;
     3057                   for (int j = 0; j < numberNonZero; j++) {
     3058                     int iSequence = index[j];
     3059                     double value = reducedCost[iSequence];
     3060                     double value2 = array[j];
     3061                     array[j] = 0.0;
     3062                     value -= value2;
     3063                     reducedCost[iSequence] = value;
     3064                     // simplify status
     3065                     ClpSimplex::Status status = model->getStatus(iSequence);
     3066         
     3067                     switch(status) {
     3068                       
     3069                     case ClpSimplex::basic:
     3070                     case ClpSimplex::isFixed:
     3071                       break;
     3072                     case ClpSimplex::isFree:
     3073                     case ClpSimplex::superBasic:
     3074                       if (fabs(value) > FREE_ACCEPT * tolerance) {
     3075                         // we are going to bias towards free (but only if reasonable)
     3076                         value *= FREE_BIAS;
     3077                         value *= value;
     3078                         // store square in list
     3079                         if (infeas[iSequence]) {
     3080                           infeas[iSequence] = value; // already there
     3081                         } else {
     3082                           array[number]=value;
     3083                           index[number++]=iSequence;
     3084                         }
     3085                       } else {
     3086                         array[number]=0.0;
     3087                         index[number++]=iSequence;
     3088                       }
     3089                       break;
     3090                     case ClpSimplex::atUpperBound:
     3091                       if (value > tolerance) {
     3092                         value *= value;
     3093                         // store square in list
     3094                         if (infeas[iSequence]) {
     3095                           infeas[iSequence] = value; // already there
     3096                         } else {
     3097                           array[number]=value;
     3098                           index[number++]=iSequence;
     3099                         }
     3100                       } else {
     3101                         array[number]=0.0;
     3102                         index[number++]=iSequence;
     3103                       }
     3104                       break;
     3105                     case ClpSimplex::atLowerBound:
     3106                       if (value < -tolerance) {
     3107                         value *= value;
     3108                         // store square in list
     3109                         if (infeas[iSequence]) {
     3110                           infeas[iSequence] = value; // already there
     3111                         } else {
     3112                           array[number]=value;
     3113                           index[number++]=iSequence;
     3114                         }
     3115                       } else {
     3116                         array[number]=0.0;
     3117                         index[number++]=iSequence;
     3118                       }
     3119                     }
     3120                   }
     3121                   dj1->setNumElements(number);
     3122                 }
    22233123               } else {
    22243124                    // use special column copy
    22253125                    // reset back
    2226                     if (killDjs)
    2227                          scaleFactor = 0.0;
    2228                     columnCopy_->transposeTimes2(model, pi, dj1, piWeight, referenceIn, devex,
     3126                    assert (!killDjs);
     3127                    //if (killDjs)
     3128                    //   scaleFactor = 0.0;
     3129                    if (infeas)
     3130                      returnCode=1;
     3131                    columnCopy_->transposeTimes2(model, pi, dj1, piWeight,
     3132                                                 infeas, reducedCost,
     3133                                                 referenceIn, devex,
    22293134                                                 reference, weights, scaleFactor);
    22303135                    numberNonZero = dj1->getNumElements();
     
    23493254     if (packed)
    23503255          dj1->setPackedMode(true);
     3256     return returnCode;
    23513257}
    23523258// Updates second array for steepest and does devex weights
     
    32534159// Creates scales for column copy (rowCopy in model may be modified)
    32544160int
    3255 ClpPackedMatrix::scale(ClpModel * model, const ClpSimplex * /*baseModel*/) const
     4161ClpPackedMatrix::scale(ClpModel * model, ClpSimplex * simplex) const
    32564162{
    32574163     //const ClpSimplex * baseModel=NULL;
     
    33924298               << smallest << largest
    33934299               << CoinMessageEol;
     4300     if (smallest *1.0e12 < largest && simplex) {
     4301       // increase tolerances
     4302       simplex->setCurrentDualTolerance(CoinMax(simplex->currentDualTolerance(),5.0e-7));
     4303       simplex->setCurrentPrimalTolerance(CoinMax(simplex->currentPrimalTolerance(),5.0e-7));
     4304     }
    33944305     if (smallest >= 0.5 && largest <= 2.0 && !deletedElements) {
    33954306          // don't bother scaling
     
    49575868{
    49585869     delete columnCopy_;
     5870     if ((model->moreSpecialOptions()&16777216)!=0) {
     5871       flags_ |= 16;
     5872       // go to exact devex (unless full steepest)
     5873       ClpPrimalColumnSteepest * pricing =
     5874         dynamic_cast<ClpPrimalColumnSteepest *>(model->primalColumnPivot());
     5875       if (pricing && pricing->mode()>1)
     5876         pricing->setMode(0);
     5877     }
    49595878     if ((flags_ & 16) != 0) {
    49605879          columnCopy_ = new ClpPackedMatrix3(model, matrix_);
     
    49785897     if (columnCopy_) {
    49795898          if (sequenceIn != -999) {
    4980                if (sequenceIn != sequenceOut) {
    4981                     if (sequenceIn < numberActiveColumns_)
    4982                          columnCopy_->swapOne(model, this, sequenceIn);
    4983                     if (sequenceOut < numberActiveColumns_)
    4984                          columnCopy_->swapOne(model, this, sequenceOut);
    4985                }
     5899            columnCopy_->swapOne(model, this, sequenceIn);
     5900            if (sequenceIn != sequenceOut)
     5901              columnCopy_->swapOne(model, this, sequenceOut);
    49865902          } else {
    49875903               // do all
     
    58696785     }
    58706786}
     6787#define ALIGNMENT 32
     6788static void * clp_align (void * memory)
     6789{
     6790  if (sizeof(int)==sizeof(void *)&&ALIGNMENT) {
     6791    CoinInt64 k = reinterpret_cast<CoinInt64> (memory);
     6792    if ((k&(ALIGNMENT-1))!=0) {
     6793      k &= ~(ALIGNMENT-1);
     6794      k += ALIGNMENT;
     6795      memory = reinterpret_cast<void *> (k);
     6796    }
     6797    return memory;
     6798  } else {
     6799    CoinInt64 k = reinterpret_cast<CoinInt64> (memory);
     6800    if ((k&(ALIGNMENT-1))!=0) {
     6801      k &= ~(ALIGNMENT-1);
     6802      k += ALIGNMENT;
     6803      memory = reinterpret_cast<void *> (k);
     6804    }
     6805    return memory;
     6806  }
     6807}
    58716808/* Default constructor. */
    58726809ClpPackedMatrix3::ClpPackedMatrix3()
    5873      : numberBlocks_(0),
    5874        numberColumns_(0),
    5875        column_(NULL),
    5876        start_(NULL),
    5877        row_(NULL),
    5878        element_(NULL),
    5879        block_(NULL)
     6810  : numberBlocks_(0),
     6811  numberColumns_(0),
     6812  numberColumnsWithGaps_(0),
     6813#if ABOCA_LITE
     6814  numberChunks_(0),
     6815#endif
     6816  numberElements_(0),
     6817  maxBlockSize_(0),
     6818  column_(NULL),
     6819  start_(NULL),
     6820  row_(NULL),
     6821  element_(NULL),
     6822  temporary_(NULL),
     6823  block_(NULL)
    58806824{
    58816825}
    58826826/* Constructor from copy. */
    58836827ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model, const CoinPackedMatrix * columnCopy)
    5884      : numberBlocks_(0),
    5885        numberColumns_(0),
    5886        column_(NULL),
    5887        start_(NULL),
    5888        row_(NULL),
    5889        element_(NULL),
    5890        block_(NULL)
    5891 {
     6828  : numberBlocks_(0),
     6829  numberColumns_(0),
     6830  numberColumnsWithGaps_(0),
     6831#if ABOCA_LITE
     6832  numberChunks_(0),
     6833#endif
     6834  numberElements_(0),
     6835  maxBlockSize_(0),
     6836  column_(NULL),
     6837  start_(NULL),
     6838  row_(NULL),
     6839  element_(NULL),
     6840  temporary_(NULL),
     6841  block_(NULL)
     6842{
     6843#ifndef COIN_AVX2
     6844#define COIN_AVX2 4
     6845#else
     6846#define HASWELL
     6847#endif
     6848#if COIN_AVX2 == 1
     6849#define COIN_AVX2_SHIFT 0
     6850#elif COIN_AVX2 == 2
     6851#define COIN_AVX2_SHIFT 1
     6852#elif COIN_AVX2 == 4
     6853#define COIN_AVX2_SHIFT 2
     6854#elif COIN_AVX2 == 8
     6855  error at present;
     6856#else
     6857  error;
     6858#endif
     6859#define COIN_ALIGN 8*COIN_AVX2 // later
     6860#define COIN_ALIGN_DOUBLE COIN_AVX2
     6861#define COIN_AVX2_CHUNK 128 // do this many at a time
    58926862#define MINBLOCK 6
    58936863#define MAXBLOCK 100
    58946864#define MAXUNROLL 10
    5895      numberColumns_ = model->getNumCols();
    5896      int numberColumns = columnCopy->getNumCols();
    5897      assert (numberColumns_ >= numberColumns);
    5898      int numberRows = columnCopy->getNumRows();
    5899      int * counts = new int[numberRows+1];
    5900      CoinZeroN(counts, numberRows + 1);
    5901      CoinBigIndex nels = 0;
    5902      CoinBigIndex nZeroEl = 0;
    5903      int iColumn;
    5904      // get matrix data pointers
    5905      const int * row = columnCopy->getIndices();
    5906      const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
    5907      const int * columnLength = columnCopy->getVectorLengths();
    5908      const double * elementByColumn = columnCopy->getElements();
    5909      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    5910           CoinBigIndex start = columnStart[iColumn];
    5911           int n = columnLength[iColumn];
    5912           CoinBigIndex end = start + n;
    5913           int kZero = 0;
    5914           for (CoinBigIndex j = start; j < end; j++) {
    5915                if(!elementByColumn[j])
    5916                     kZero++;
    5917           }
    5918           n -= kZero;
    5919           nZeroEl += kZero;
    5920           nels += n;
    5921           counts[n]++;
    5922      }
    5923      counts[0] += numberColumns_ - numberColumns;
    5924      int nZeroColumns = counts[0];
    5925      counts[0] = -1;
    5926      numberColumns_ -= nZeroColumns;
    5927      column_ = new int [2*numberColumns_+nZeroColumns];
    5928      int * lookup = column_ + numberColumns_;
    5929      row_ = new int[nels];
    5930      element_ = new double[nels];
    5931      int nOdd = 0;
    5932      CoinBigIndex nInOdd = 0;
    5933      int i;
    5934      for (i = 1; i <= numberRows; i++) {
    5935           int n = counts[i];
    5936           if (n) {
    5937                if (n < MINBLOCK || i > MAXBLOCK) {
    5938                     nOdd += n;
    5939                     counts[i] = -1;
    5940                     nInOdd += n * i;
    5941                } else {
    5942                     numberBlocks_++;
    5943                }
    5944           } else {
    5945                counts[i] = -1;
    5946           }
    5947      }
    5948      start_ = new CoinBigIndex [nOdd+1];
    5949      // even if no blocks do a dummy one
    5950      numberBlocks_ = CoinMax(numberBlocks_, 1);
    5951      block_ = new blockStruct [numberBlocks_];
    5952      memset(block_, 0, numberBlocks_ * sizeof(blockStruct));
    5953      // Fill in what we can
    5954      int nTotal = nOdd;
    5955      // in case no blocks
    5956      block_->startIndices_ = nTotal;
    5957      nels = nInOdd;
    5958      int nBlock = 0;
    5959      for (i = 0; i <= CoinMin(MAXBLOCK, numberRows); i++) {
    5960           if (counts[i] > 0) {
    5961                blockStruct * block = block_ + nBlock;
    5962                int n = counts[i];
    5963                counts[i] = nBlock; // backward pointer
    5964                nBlock++;
    5965                block->startIndices_ = nTotal;
    5966                block->startElements_ = nels;
    5967                block->numberElements_ = i;
    5968                // up counts
    5969                nTotal += n;
    5970                nels += n * i;
    5971           }
    5972      }
    5973      for (iColumn = numberColumns; iColumn < numberColumns_; iColumn++)
    5974           lookup[iColumn] = -1;
    5975      // fill
    5976      start_[0] = 0;
    5977      nOdd = 0;
    5978      nInOdd = 0;
    5979      const double * columnScale = model->columnScale();
    5980      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    5981           CoinBigIndex start = columnStart[iColumn];
    5982           int n = columnLength[iColumn];
    5983           CoinBigIndex end = start + n;
    5984           int kZero = 0;
    5985           for (CoinBigIndex j = start; j < end; j++) {
    5986                if(!elementByColumn[j])
    5987                     kZero++;
    5988           }
    5989           n -= kZero;
    5990           if (n) {
    5991                int iBlock = counts[n];
    5992                if (iBlock >= 0) {
    5993                     blockStruct * block = block_ + iBlock;
    5994                     int k = block->numberInBlock_;
    5995                     block->numberInBlock_ ++;
    5996                     column_[block->startIndices_+k] = iColumn;
    5997                     lookup[iColumn] = k;
    5998                     CoinBigIndex put = block->startElements_ + k * n;
    5999                     for (CoinBigIndex j = start; j < end; j++) {
    6000                          double value = elementByColumn[j];
    6001                          if(value) {
    6002                               if (columnScale)
    6003                                    value *= columnScale[iColumn];
    6004                               element_[put] = value;
    6005                               row_[put++] = row[j];
    6006                          }
    6007                     }
    6008                } else {
    6009                     // odd ones
    6010                     for (CoinBigIndex j = start; j < end; j++) {
    6011                          double value = elementByColumn[j];
    6012                          if(value) {
    6013                               if (columnScale)
    6014                                    value *= columnScale[iColumn];
    6015                               element_[nInOdd] = value;
    6016                               row_[nInOdd++] = row[j];
    6017                          }
    6018                     }
    6019                     column_[nOdd] = iColumn;
    6020                     lookup[iColumn] = -1;
    6021                     nOdd++;
    6022                     start_[nOdd] = nInOdd;
    6023                }
    6024           } else {
    6025                // zero column
    6026                lookup[iColumn] = -1;
    6027           }
    6028      }
    6029      delete [] counts;
     6865#define roundUp(xx) ((xx+COIN_ALIGN_DOUBLE-1)&(~(COIN_ALIGN_DOUBLE-1)))
     6866#define roundDown(xx) (xx&(~(COIN_ALIGN_DOUBLE-1)))
     6867  numberColumns_ = model->getNumCols();
     6868  int numberColumns = columnCopy->getNumCols();
     6869  assert (numberColumns_ >= numberColumns);
     6870  int numberRows = columnCopy->getNumRows();
     6871  int * counts = new int[numberRows+1];
     6872  CoinZeroN(counts, numberRows + 1);
     6873  CoinBigIndex nels = 0;
     6874  int iColumn;
     6875  // get matrix data pointers
     6876  const int * row = columnCopy->getIndices();
     6877  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     6878  const int * columnLength = columnCopy->getVectorLengths();
     6879  const double * elementByColumn = columnCopy->getElements();
     6880  unsigned char * status = model->statusArray();
     6881  const double * lower = model->columnLower();
     6882  const double * upper = model->columnUpper();
     6883  int nFree=0; // or superbasic
     6884  int nFreeEls=0;
     6885  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     6886    CoinBigIndex start = columnStart[iColumn];
     6887    int n = columnLength[iColumn];
     6888    CoinBigIndex end = start + n;
     6889    int kZero = 0;
     6890    for (CoinBigIndex j = start; j < end; j++) {
     6891      if(!elementByColumn[j])
     6892        kZero++;
     6893    }
     6894    n -= kZero;
     6895    nels += n;
     6896    if ((lower[iColumn]==-COIN_DBL_MAX&&upper[iColumn]==COIN_DBL_MAX)||
     6897        (status[iColumn]&3) == 0) {
     6898      nFree++;
     6899      nFreeEls += n;
     6900      n=0;
     6901      if ((status[iColumn]&3) != 0) {
     6902        status[iColumn] &= ~7;
     6903        status[iColumn] |= 4;
     6904      }
     6905    }
     6906    counts[n]++;
     6907  }
     6908  counts[0] += numberColumns_ - numberColumns;
     6909  int nZeroColumns = counts[0]; // also free
     6910  counts[0] = -1;
     6911  int nOdd = nZeroColumns;
     6912  CoinBigIndex nInOdd = nFreeEls;
     6913  maxBlockSize_ = 0;
     6914  int i;
     6915  for (i = 1; i <= numberRows; i++) {
     6916    int n = counts[i];
     6917    if (n) {
     6918      if (n < MINBLOCK || i > MAXBLOCK) {
     6919        nOdd += n;
     6920        counts[i] = -1;
     6921        nInOdd += n * i;
     6922      } else {
     6923        numberBlocks_++;
     6924        maxBlockSize_ = CoinMax(maxBlockSize_,n);;
     6925      }
     6926    } else {
     6927      counts[i] = -1;
     6928    }
     6929  }
     6930  // align on boundaries
     6931  nels = roundUp(nInOdd);
     6932  numberColumnsWithGaps_ = nOdd;
     6933  for (int i = 0; i <= CoinMin(MAXBLOCK, numberRows); i++) {
     6934    if (counts[i] > 0) {
     6935      int n = roundUp(counts[i]);
     6936      nels += n*i;
     6937      numberColumnsWithGaps_ += n;
     6938    }
     6939  }
     6940  row_ = new int[nels+15];
     6941  element_ = new double[nels+31];
     6942  start_ = new CoinBigIndex [nOdd+1];
     6943  // also allow for rows!
     6944  int numberColumnsWithGaps = roundUp(numberColumnsWithGaps_);
     6945  numberColumnsWithGaps_ = roundUp(numberColumnsWithGaps+numberRows);
     6946  column_ = new int [2*numberColumnsWithGaps_];
     6947  memset(row_,0,nels*sizeof(int));
     6948  memset(element_,0,nels*sizeof(double));
     6949  int * lookup = column_ + numberColumnsWithGaps_;
     6950  for (int i = 0; i < numberColumnsWithGaps; i++) {
     6951    column_[i] = -1;
     6952    lookup[i] = -1;
     6953  }
     6954  for (int i = 0; i < numberRows; i++) {
     6955    column_[i+numberColumnsWithGaps] = i+numberColumns;
     6956    lookup[i+numberColumns] = i;
     6957  }
     6958  for (int i = numberRows+numberColumnsWithGaps;
     6959       i < numberColumnsWithGaps_; i++) {
     6960    column_[i] = -1;
     6961    lookup[i] = -1;
     6962  }
     6963  // even if no blocks do a dummy one
     6964  numberBlocks_ = CoinMax(numberBlocks_, 1);
     6965  block_ = new blockStruct [numberBlocks_+1]; // +1 for slacks
     6966  memset(block_, 0, (numberBlocks_+1) * sizeof(blockStruct));
     6967  // Fill in what we can
     6968  int nTotal = nOdd;
     6969  // in case no blocks
     6970  block_->startIndices_ = nTotal;
     6971  // adjust nels so start on 8 double boundary
     6972  double * elements2 = reinterpret_cast<double *>(clp_align(element_+nInOdd));
     6973  nels = elements2-element_;
     6974  int nBlock = 0;
     6975  for (i = 0; i <= CoinMin(MAXBLOCK, numberRows); i++) {
     6976    if (counts[i] > 0) {
     6977      blockStruct * block = block_ + nBlock;
     6978      int n = roundUp(counts[i]);
     6979      counts[i] = nBlock; // backward pointer
     6980      nBlock++;
     6981      block->startIndices_ = nTotal;
     6982      block->startElements_ = nels;
     6983      block->numberElements_ = i;
     6984      // up counts
     6985      nTotal += n;
     6986      nels += n * i;
     6987    }
     6988  }
     6989  numberElements_ = nels;
     6990  nBlock = CoinMax(nBlock, 1);
     6991  // slacks
     6992  block_[nBlock].numberElements_ = 0;
     6993  block_[nBlock].numberInBlock_ = numberRows;
     6994  block_[nBlock].startIndices_ = numberColumnsWithGaps;
     6995  // fill
     6996  start_[0] = 0;
     6997  nOdd = 0;
     6998  nInOdd = 0;
     6999  const double * columnScale = model->columnScale();
     7000  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     7001    CoinBigIndex start = columnStart[iColumn];
     7002    int n = columnLength[iColumn];
     7003    CoinBigIndex end = start + n;
     7004    int kZero = 0;
     7005    for (CoinBigIndex j = start; j < end; j++) {
     7006      if(!elementByColumn[j])
     7007        kZero++;
     7008    }
     7009    n -= kZero;
     7010    if ((status[iColumn]&3) == 0)
     7011      n=0;
     7012    {
     7013      int iBlock = counts[n];
     7014      if (iBlock >= 0) {
     7015        blockStruct * block = block_ + iBlock;
     7016        int k = block->numberInBlock_;
     7017        block->numberInBlock_ ++;
     7018        column_[block->startIndices_+k] = iColumn;
     7019        lookup[iColumn] = k;
     7020        CoinBigIndex put = block->startElements_
     7021          + roundDown(k) * n + (k&(COIN_AVX2-1));
     7022        for (CoinBigIndex j = start; j < end; j++) {
     7023          double value = elementByColumn[j];
     7024          if(value) {
     7025            if (columnScale)
     7026              value *= columnScale[iColumn];
     7027            element_[put] = value;
     7028            row_[put] = row[j];
     7029            put += COIN_AVX2;
     7030          }
     7031        }
     7032      } else {
     7033        // odd ones
     7034        for (CoinBigIndex j = start; j < end; j++) {
     7035          double value = elementByColumn[j];
     7036          if(value) {
     7037            if (columnScale)
     7038              value *= columnScale[iColumn];
     7039            element_[nInOdd] = value;
     7040            row_[nInOdd++] = row[j];
     7041          }
     7042        }
     7043        column_[nOdd] = iColumn;
     7044        lookup[iColumn] = -1;
     7045        nOdd++;
     7046        start_[nOdd] = nInOdd;
     7047      }
     7048    }
     7049  }
     7050  // This much temporary space
     7051#if ABOCA_LITE
     7052  temporary_ = new CoinDoubleArrayWithLength(2*COIN_AVX2_CHUNK*2*ABOCA_LITE,
     7053                                             -2-COIN_ALIGN_DOUBLE);
     7054#else
     7055  temporary_ = new CoinDoubleArrayWithLength(2*COIN_AVX2_CHUNK,
     7056                                             -2-COIN_ALIGN_DOUBLE);
     7057#endif
     7058#if ABOCA_LITE
     7059  // balance work - maybe with 2*cpu chunks
     7060#define COLUMN_WEIGHT 5
     7061#define BLOCK_MULT 2
     7062  // If do +1's then need two element weights
     7063  double totalWork=0.0;
     7064  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     7065    totalWork += (COLUMN_WEIGHT+block_[iBlock].numberElements_)*
     7066      block_[iBlock].numberInBlock_;
     7067    if (model->logLevel()>2)
     7068    printf("block %d ncol %d nel %d total %g this %d\n",
     7069           iBlock,block_[iBlock].numberInBlock_,
     7070           block_[iBlock].numberElements_,totalWork,
     7071           (COLUMN_WEIGHT+block_[iBlock].numberElements_)*
     7072           block_[iBlock].numberInBlock_);
     7073  }
     7074  double eachWork = totalWork/(BLOCK_MULT*ABOCA_LITE);
     7075  double thisWork = 0.0;
     7076  numberChunks_=0;
     7077  endChunk_[0]=0;
     7078  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     7079    thisWork += (COLUMN_WEIGHT+block_[iBlock].numberElements_)*
     7080      block_[iBlock].numberInBlock_;
     7081    if (thisWork>=eachWork) {
     7082      totalWork -= thisWork;
     7083      thisWork=0.0;
     7084      eachWork = totalWork/(BLOCK_MULT*ABOCA_LITE-numberChunks_-1);
     7085      numberChunks_++;
     7086      endChunk_[numberChunks_]=iBlock+1;
     7087      if (numberChunks_==BLOCK_MULT*ABOCA_LITE)
     7088        break;
     7089    }
     7090  }
     7091  if (numberChunks_==BLOCK_MULT*ABOCA_LITE)
     7092    endChunk_[numberChunks_]=numberBlocks_;
     7093  else if (endChunk_[numberChunks_]<numberBlocks_)
     7094    endChunk_[++numberChunks_]=numberBlocks_;
     7095  //printf("%d chunks for %d blocks\n",numberChunks_,numberBlocks_);
     7096  assert(numberChunks_<=2*ABOCA_LITE);
     7097#endif
     7098  delete [] counts;
    60307099}
    60317100/* Destructor */
    60327101ClpPackedMatrix3::~ClpPackedMatrix3()
    60337102{
    6034      delete [] column_;
    6035      delete [] start_;
    6036      delete [] row_;
    6037      delete [] element_;
    6038      delete [] block_;
     7103  delete [] column_;
     7104  delete [] start_;
     7105  delete [] row_;
     7106  delete [] element_;
     7107  delete temporary_;
     7108  delete [] block_;
    60397109}
    60407110/* The copy constructor. */
    60417111ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
    6042      : numberBlocks_(rhs.numberBlocks_),
    6043        numberColumns_(rhs.numberColumns_),
    6044        column_(NULL),
    6045        start_(NULL),
    6046        row_(NULL),
    6047        element_(NULL),
    6048        block_(NULL)
    6049 {
    6050      if (rhs.numberBlocks_) {
    6051           block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
    6052           column_ = CoinCopyOfArray(rhs.column_, 2 * numberColumns_);
    6053           int numberOdd = block_->startIndices_;
    6054           start_ = CoinCopyOfArray(rhs.start_, numberOdd + 1);
    6055           blockStruct * lastBlock = block_ + (numberBlocks_ - 1);
    6056           CoinBigIndex numberElements = lastBlock->startElements_ +
    6057                                         lastBlock->numberInBlock_ * lastBlock->numberElements_;
    6058           row_ = CoinCopyOfArray(rhs.row_, numberElements);
    6059           element_ = CoinCopyOfArray(rhs.element_, numberElements);
    6060      }
     7112  : numberBlocks_(rhs.numberBlocks_),
     7113  numberColumns_(rhs.numberColumns_),
     7114  numberColumnsWithGaps_(rhs.numberColumnsWithGaps_),
     7115#if ABOCA_LITE
     7116  numberChunks_(rhs.numberChunks_),
     7117#endif
     7118  numberElements_(rhs.numberElements_),
     7119  maxBlockSize_(rhs.maxBlockSize_),
     7120  column_(NULL),
     7121  start_(NULL),
     7122  row_(NULL),
     7123  element_(NULL),
     7124  temporary_(NULL),
     7125  block_(NULL)
     7126{
     7127  if (rhs.numberBlocks_) {
     7128    block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
     7129    column_ = CoinCopyOfArray(rhs.column_, 2 * numberColumnsWithGaps_);
     7130    int numberOdd = block_->startIndices_;
     7131    start_ = CoinCopyOfArray(rhs.start_, numberOdd + 1);
     7132    row_ = CoinCopyOfArray(rhs.row_, numberElements_);
     7133    element_ = CoinCopyOfArray(rhs.element_, numberElements_+8);
     7134    // This much temporary space
     7135#if ABOCA_LITE
     7136    temporary_ = new CoinDoubleArrayWithLength(2*COIN_AVX2_CHUNK*2*ABOCA_LITE,
     7137                                               -2-COIN_ALIGN_DOUBLE);
     7138#else
     7139    temporary_ = new CoinDoubleArrayWithLength(2*COIN_AVX2_CHUNK,
     7140                                               -2-COIN_ALIGN_DOUBLE);
     7141#endif
     7142#if ABOCA_LITE
     7143    memcpy(endChunk_,rhs.endChunk_,sizeof(endChunk_));
     7144#endif
     7145  }
    60617146}
    60627147ClpPackedMatrix3&
    6063 ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
    6064 {
    6065      if (this != &rhs) {
    6066           delete [] column_;
    6067           delete [] start_;
    6068           delete [] row_;
    6069           delete [] element_;
    6070           delete [] block_;
    6071           numberBlocks_ = rhs.numberBlocks_;
    6072           numberColumns_ = rhs.numberColumns_;
    6073           if (rhs.numberBlocks_) {
    6074                block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
    6075                column_ = CoinCopyOfArray(rhs.column_, 2 * numberColumns_);
    6076                int numberOdd = block_->startIndices_;
    6077                start_ = CoinCopyOfArray(rhs.start_, numberOdd + 1);
    6078                blockStruct * lastBlock = block_ + (numberBlocks_ - 1);
    6079                CoinBigIndex numberElements = lastBlock->startElements_ +
    6080                                              lastBlock->numberInBlock_ * lastBlock->numberElements_;
    6081                row_ = CoinCopyOfArray(rhs.row_, numberElements);
    6082                element_ = CoinCopyOfArray(rhs.element_, numberElements);
    6083           } else {
    6084                column_ = NULL;
    6085                start_ = NULL;
    6086                row_ = NULL;
    6087                element_ = NULL;
    6088                block_ = NULL;
    6089           }
    6090      }
    6091      return *this;
     7148   ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
     7149{
     7150  if (this != &rhs) {
     7151    delete [] column_;
     7152    delete [] start_;
     7153    delete [] row_;
     7154    delete [] element_;
     7155    delete temporary_;
     7156    temporary_=NULL;
     7157    delete [] block_;
     7158    numberBlocks_ = rhs.numberBlocks_;
     7159    numberColumns_ = rhs.numberColumns_;
     7160    numberColumnsWithGaps_ = rhs.numberColumnsWithGaps_;
     7161#if ABOCA_LITE
     7162    numberChunks_ = rhs.numberChunks_;
     7163    memcpy(endChunk_,rhs.endChunk_,sizeof(endChunk_));
     7164#endif
     7165    numberElements_ = rhs.numberElements_;
     7166    maxBlockSize_ = rhs.maxBlockSize_;
     7167    if (rhs.numberBlocks_) {
     7168      block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
     7169      column_ = CoinCopyOfArray(rhs.column_, 2 * numberColumnsWithGaps_);
     7170      int numberOdd = block_->startIndices_;
     7171      start_ = CoinCopyOfArray(rhs.start_, numberOdd + 1);
     7172      row_ = CoinCopyOfArray(rhs.row_, numberElements_);
     7173      element_ = CoinCopyOfArray(rhs.element_, numberElements_+8);
     7174      // This much temporary space
     7175#if ABOCA_LITE
     7176      temporary_ = new CoinDoubleArrayWithLength(2*COIN_AVX2_CHUNK*2*ABOCA_LITE,
     7177                                                 -2-COIN_ALIGN_DOUBLE);
     7178#else
     7179      temporary_ = new CoinDoubleArrayWithLength(2*COIN_AVX2_CHUNK,
     7180                                                 -2-COIN_ALIGN_DOUBLE);
     7181#endif
     7182    } else {
     7183      column_ = NULL;
     7184      start_ = NULL;
     7185      row_ = NULL;
     7186      element_ = NULL;
     7187      block_ = NULL;
     7188    }
     7189  }
     7190  return *this;
    60927191}
    60937192/* Sort blocks */
     
    60957194ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
    60967195{
    6097      int * lookup = column_ + numberColumns_;
    6098      for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
    6099           blockStruct * block = block_ + iBlock;
    6100           int numberInBlock = block->numberInBlock_;
    6101           int nel = block->numberElements_;
    6102           int * row = row_ + block->startElements_;
    6103           double * element = element_ + block->startElements_;
    6104           int * column = column_ + block->startIndices_;
    6105           int lastPrice = 0;
    6106           int firstNotPrice = numberInBlock - 1;
    6107           while (lastPrice <= firstNotPrice) {
    6108                // find first basic or fixed
    6109                int iColumn = numberInBlock;
    6110                for (; lastPrice <= firstNotPrice; lastPrice++) {
    6111                     iColumn = column[lastPrice];
    6112                     if (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
    6113                               model->getColumnStatus(iColumn) == ClpSimplex::isFixed)
    6114                          break;
    6115                }
    6116                // find last non basic or fixed
    6117                int jColumn = -1;
    6118                for (; firstNotPrice > lastPrice; firstNotPrice--) {
    6119                     jColumn = column[firstNotPrice];
    6120                     if (model->getColumnStatus(jColumn) != ClpSimplex::basic &&
    6121                               model->getColumnStatus(jColumn) != ClpSimplex::isFixed)
    6122                          break;
    6123                }
    6124                if (firstNotPrice > lastPrice) {
    6125                     assert (column[lastPrice] == iColumn);
    6126                     assert (column[firstNotPrice] == jColumn);
    6127                     // need to swap
    6128                     column[firstNotPrice] = iColumn;
    6129                     lookup[iColumn] = firstNotPrice;
    6130                     column[lastPrice] = jColumn;
    6131                     lookup[jColumn] = lastPrice;
    6132                     double * elementA = element + lastPrice * nel;
    6133                     int * rowA = row + lastPrice * nel;
    6134                     double * elementB = element + firstNotPrice * nel;
    6135                     int * rowB = row + firstNotPrice * nel;
    6136                     for (int i = 0; i < nel; i++) {
    6137                          int temp = rowA[i];
    6138                          double tempE = elementA[i];
    6139                          rowA[i] = rowB[i];
    6140                          elementA[i] = elementB[i];
    6141                          rowB[i] = temp;
    6142                          elementB[i] = tempE;
    6143                     }
    6144                     firstNotPrice--;
    6145                     lastPrice++;
    6146                } else if (lastPrice == firstNotPrice) {
    6147                     // make sure correct side
    6148                     iColumn = column[lastPrice];
    6149                     if (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
    6150                               model->getColumnStatus(iColumn) != ClpSimplex::isFixed)
    6151                          lastPrice++;
    6152                     break;
    6153                }
    6154           }
    6155           block->numberPrice_ = lastPrice;
     7196  int * lookup = column_ + numberColumnsWithGaps_;
     7197  for (int iBlock = 0; iBlock < numberBlocks_+1; iBlock++) {
     7198    blockStruct * block = block_ + iBlock;
     7199    int numberInBlock = block->numberInBlock_;
     7200    int nel = block->numberElements_;
     7201    int * row = row_ + block->startElements_;
     7202    double * element = element_ + block->startElements_;
     7203    int * column = column_ + block->startIndices_;
     7204    int lastPrice = 0;
     7205    int firstNotPrice = numberInBlock - 1;
     7206    while (lastPrice <= firstNotPrice) {
     7207      // find first basic or fixed
     7208      int iColumn = numberInBlock;
     7209      for (; lastPrice <= firstNotPrice; lastPrice++) {
     7210        iColumn = column[lastPrice];
     7211        if (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
     7212            model->getColumnStatus(iColumn) == ClpSimplex::isFixed)
     7213          break;
     7214      }
     7215      // find last non basic or fixed
     7216      int jColumn = -1;
     7217      for (; firstNotPrice > lastPrice; firstNotPrice--) {
     7218        jColumn = column[firstNotPrice];
     7219        if (model->getColumnStatus(jColumn) != ClpSimplex::basic &&
     7220            model->getColumnStatus(jColumn) != ClpSimplex::isFixed)
     7221          break;
     7222      }
     7223      if (firstNotPrice > lastPrice) {
     7224        assert (column[lastPrice] == iColumn);
     7225        assert (column[firstNotPrice] == jColumn);
     7226        // need to swap
     7227        column[firstNotPrice] = iColumn;
     7228        lookup[iColumn] = firstNotPrice;
     7229        column[lastPrice] = jColumn;
     7230        lookup[jColumn] = lastPrice;
     7231        int startBit = roundDown(lastPrice);
     7232        CoinBigIndex offset =startBit*nel + (lastPrice-startBit);
     7233        double * elementA = element + offset;
     7234        int * rowA = row + offset;
     7235        startBit = roundDown(firstNotPrice);
     7236        offset =startBit*nel + (firstNotPrice-startBit);
     7237        double * elementB = element + offset;
     7238        int * rowB = row + offset;
     7239        for (int i = 0; i < nel*COIN_AVX2; i+=COIN_AVX2) {
     7240          int temp = rowA[i];
     7241          double tempE = elementA[i];
     7242          rowA[i] = rowB[i];
     7243          elementA[i] = elementB[i];
     7244          rowB[i] = temp;
     7245          elementB[i] = tempE;
     7246        }
     7247        firstNotPrice--;
     7248        lastPrice++;
     7249      } else if (lastPrice == firstNotPrice) {
     7250        // make sure correct side
     7251        iColumn = column[lastPrice];
     7252        if (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
     7253            model->getColumnStatus(iColumn) != ClpSimplex::isFixed)
     7254          lastPrice++;
     7255        break;
     7256      }
     7257    }
     7258    block->firstBasic_ = lastPrice;
     7259    // reuse ints even though names wrong
     7260    // lower
     7261    firstNotPrice = lastPrice - 1;
     7262    lastPrice = 0;
     7263    while (lastPrice <= firstNotPrice) {
     7264      // find first lower
     7265      int iColumn = numberInBlock;
     7266      for (; lastPrice <= firstNotPrice; lastPrice++) {
     7267        iColumn = column[lastPrice];
     7268        if (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound)
     7269          break;
     7270      }
     7271      // find last lower
     7272      int jColumn = -1;
     7273      for (; firstNotPrice > lastPrice; firstNotPrice--) {
     7274        jColumn = column[firstNotPrice];
     7275        if (model->getColumnStatus(jColumn) != ClpSimplex::atLowerBound)
     7276          break;
     7277      }
     7278      if (firstNotPrice > lastPrice) {
     7279        assert (column[lastPrice] == iColumn);
     7280        assert (column[firstNotPrice] == jColumn);
     7281        // need to swap
     7282        column[firstNotPrice] = iColumn;
     7283        lookup[iColumn] = firstNotPrice;
     7284        column[lastPrice] = jColumn;
     7285        lookup[jColumn] = lastPrice;
     7286        int startBit = roundDown(lastPrice);
     7287        CoinBigIndex offset =startBit*nel + (lastPrice-startBit);
     7288        double * elementA = element + offset;
     7289        int * rowA = row + offset;
     7290        startBit = roundDown(firstNotPrice);
     7291        offset =startBit*nel + (firstNotPrice-startBit);
     7292        double * elementB = element + offset;
     7293        int * rowB = row + offset;
     7294        for (int i = 0; i < nel*COIN_AVX2; i+=COIN_AVX2) {
     7295          int temp = rowA[i];
     7296          double tempE = elementA[i];
     7297          rowA[i] = rowB[i];
     7298          elementA[i] = elementB[i];
     7299          rowB[i] = temp;
     7300          elementB[i] = tempE;
     7301        }
     7302        firstNotPrice--;
     7303        lastPrice++;
     7304      } else if (lastPrice == firstNotPrice) {
     7305        // make sure correct side
     7306        iColumn = column[lastPrice];
     7307        if (model->getColumnStatus(iColumn) != ClpSimplex::atLowerBound)
     7308          lastPrice++;
     7309        break;
     7310      }
     7311    }
     7312    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;
    61567366#ifndef NDEBUG
    6157           // check
    6158           int i;
    6159           for (i = 0; i < lastPrice; i++) {
    6160                int iColumn = column[i];
    6161                assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
    6162                        model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
    6163                assert (lookup[iColumn] == i);
    6164           }
    6165           for (; i < numberInBlock; i++) {
    6166                int iColumn = column[i];
    6167                assert (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
    6168                        model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
    6169                assert (lookup[iColumn] == i);
    6170           }
    6171 #endif
    6172      }
     7367    // check
     7368    int i;
     7369    for (i = 0; i < block->firstBasic_; i++) {
     7370      int iColumn = column[i];
     7371      assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
     7372              model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
     7373      assert (lookup[iColumn] == i);
     7374      if (i<block->firstAtUpper_) {
     7375        assert (model->getColumnStatus(iColumn) == ClpSimplex::isFree ||
     7376                model->getColumnStatus(iColumn) == ClpSimplex::superBasic);
     7377      } else if (i<block->firstAtLower_) {
     7378        assert (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound);
     7379      } else {
     7380        assert (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound);
     7381      }
     7382    }
     7383    for (; i < numberInBlock; i++) {
     7384      int iColumn = column[i];
     7385      assert (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
     7386              model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
     7387      assert (lookup[iColumn] == i);
     7388    }
     7389#endif
     7390  }
     7391}
     7392// Swap one variable - when found
     7393void
     7394  ClpPackedMatrix3::swapOne(int iBlock, int kA, int kB)
     7395{
     7396  int * lookup = column_ + numberColumnsWithGaps_;
     7397  blockStruct * block = block_ + iBlock;
     7398  int nel = block->numberElements_;
     7399  int * row = row_ + block->startElements_;
     7400  double * element = element_ + block->startElements_;
     7401  int * column = column_ + block->startIndices_;
     7402  int iColumn = column[kA];
     7403  int jColumn = column[kB];
     7404  column[kA] = jColumn;
     7405  lookup[jColumn] = kA;
     7406  column[kB] = iColumn;
     7407  lookup[iColumn] = kB;
     7408  int startBit = roundDown(kA);
     7409  CoinBigIndex offset =startBit*nel + (kA-startBit);
     7410  double * elementA = element + offset;
     7411  int * rowA = row + offset;
     7412  startBit = roundDown(kB);
     7413  offset =startBit*nel + (kB-startBit);
     7414  double * elementB = element + offset;
     7415  int * rowB = row + offset;
     7416  int i;
     7417  for (i = 0; i < nel*COIN_AVX2; i+=COIN_AVX2) {
     7418    int temp = rowA[i];
     7419    double tempE = elementA[i];
     7420    rowA[i] = rowB[i];
     7421    elementA[i] = elementB[i];
     7422    rowB[i] = temp;
     7423    elementB[i] = tempE;
     7424  }
    61737425}
    61747426// Swap one variable
    61757427void
    6176 ClpPackedMatrix3::swapOne(const ClpSimplex * model, const ClpPackedMatrix * matrix,
    6177                           int iColumn)
    6178 {
    6179      int * lookup = column_ + numberColumns_;
    6180      // position in block
    6181      int kA = lookup[iColumn];
    6182      if (kA < 0)
    6183           return; // odd one
    6184      // get matrix data pointers
    6185      const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
    6186      //const int * row = columnCopy->getIndices();
    6187      const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
    6188      const int * columnLength = columnCopy->getVectorLengths();
    6189      const double * elementByColumn = columnCopy->getElements();
    6190      CoinBigIndex start = columnStart[iColumn];
    6191      int n = columnLength[iColumn];
    6192      if (matrix->zeros()) {
    6193           CoinBigIndex end = start + n;
    6194           for (CoinBigIndex j = start; j < end; j++) {
    6195                if(!elementByColumn[j])
    6196                     n--;
    6197           }
    6198      }
    6199      // find block - could do binary search
    6200      int iBlock = CoinMin(n, numberBlocks_) - 1;
    6201      while (block_[iBlock].numberElements_ != n)
    6202           iBlock--;
    6203      blockStruct * block = block_ + iBlock;
    6204      int nel = block->numberElements_;
    6205      int * row = row_ + block->startElements_;
    6206      double * element = element_ + block->startElements_;
    6207      int * column = column_ + block->startIndices_;
    6208      assert (column[kA] == iColumn);
    6209      bool moveUp = (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
    6210                     model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
    6211      int lastPrice = block->numberPrice_;
    6212      int kB;
    6213      if (moveUp) {
    6214           // May already be in correct place (e.g. fixed basic leaving basis)
    6215           if (kA >= lastPrice)
    6216                return;
    6217           kB = lastPrice - 1;
    6218           block->numberPrice_--;
    6219      } else {
    6220           assert (kA >= lastPrice);
    6221           kB = lastPrice;
    6222           block->numberPrice_++;
    6223      }
    6224      int jColumn = column[kB];
    6225      column[kA] = jColumn;
    6226      lookup[jColumn] = kA;
    6227      column[kB] = iColumn;
    6228      lookup[iColumn] = kB;
    6229      double * elementA = element + kB * nel;
    6230      int * rowA = row + kB * nel;
    6231      double * elementB = element + kA * nel;
    6232      int * rowB = row + kA * nel;
    6233      int i;
    6234      for (i = 0; i < nel; i++) {
    6235           int temp = rowA[i];
    6236           double tempE = elementA[i];
    6237           rowA[i] = rowB[i];
    6238           elementA[i] = elementB[i];
    6239           rowB[i] = temp;
    6240           elementB[i] = tempE;
    6241      }
     7428  ClpPackedMatrix3::swapOne(const ClpSimplex * model, const ClpPackedMatrix * matrix,
     7429                            int iColumn)
     7430{
     7431  int * lookup = column_ + numberColumnsWithGaps_;
     7432  // position in block
     7433  int kA = lookup[iColumn];
     7434  if (kA < 0)
     7435    return; // odd one
     7436  int iBlock = numberBlocks_;
     7437  if (iColumn<model->numberColumns()) {
     7438    // get matrix data pointers
     7439    const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
     7440    //const int * row = columnCopy->getIndices();
     7441    const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     7442    const int * columnLength = columnCopy->getVectorLengths();
     7443    const double * elementByColumn = columnCopy->getElements();
     7444    CoinBigIndex start = columnStart[iColumn];
     7445    int n = columnLength[iColumn];
     7446    if (matrix->zeros()) {
     7447      CoinBigIndex end = start + n;
     7448      for (CoinBigIndex j = start; j < end; j++) {
     7449        if(!elementByColumn[j])
     7450          n--;
     7451      }
     7452    }
     7453    // find block - could do binary search
     7454    iBlock = CoinMin(n, numberBlocks_) - 1;
     7455    while (block_[iBlock].numberElements_ != n)
     7456      iBlock--;
     7457  }
     7458  blockStruct * block = block_ + iBlock;
    62427459#ifndef NDEBUG
    6243      // check
    6244      for (i = 0; i < block->numberPrice_; i++) {
    6245           int iColumn = column[i];
    6246           if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
    6247                assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
    6248                        model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
    6249           assert (lookup[iColumn] == i);
    6250      }
    6251      int numberInBlock = block->numberInBlock_;
    6252      for (; i < numberInBlock; i++) {
    6253           int iColumn = column[i];
    6254           if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
    6255                assert (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
    6256                        model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
    6257           assert (lookup[iColumn] == i);
    6258      }
     7460  int * column = column_ + block->startIndices_;
     7461  assert (column[kA] == iColumn);
     7462#endif
     7463  int from,to;
     7464  if (kA<block->firstBasic_) {
     7465    if (kA>=block->firstAtLower_) {
     7466      from=2;
     7467    } else if (kA>=block->firstAtUpper_) {
     7468      from=1;
     7469    } else {
     7470      from=0;
     7471    }
     7472  } else {
     7473    from=3;
     7474  }
     7475  if (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
     7476      model->getColumnStatus(iColumn) == ClpSimplex::isFixed) {
     7477    to=3;
     7478  } else if (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound) {
     7479    to=2;
     7480  } else if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
     7481    to=1;
     7482  } else {
     7483    to=0;
     7484  }
     7485  int * statusCounts = (&block->firstAtUpper_)-1;
     7486  if (from<to) {
     7487    while (from<to) {
     7488      int kB=statusCounts[from+1]-1;
     7489      statusCounts[from+1]=kB;
     7490      swapOne(iBlock,kA,kB);
     7491      kA=kB;
     7492      from++;
     7493    }
     7494  } else if (from>to) {
     7495    while (from>to) {
     7496      int kB=statusCounts[from];
     7497      statusCounts[from]=kB+1;
     7498      swapOne(iBlock,kA,kB);
     7499      kA=kB;
     7500      from--;
     7501    }
     7502  }
     7503#ifndef NDEBUG
     7504  // check
     7505  int i;
     7506  for (i = 0; i < block->firstBasic_; i++) {
     7507    int iColumn = column[i];
     7508    if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
     7509      assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
     7510              model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
     7511    assert (lookup[iColumn] == i);
     7512    if (model->algorithm()>0) {
     7513      if (i<block->firstAtUpper_) {
     7514        assert (model->getColumnStatus(iColumn) == ClpSimplex::isFree ||
     7515                model->getColumnStatus(iColumn) == ClpSimplex::superBasic);
     7516      } else if (i<block->firstAtLower_) {
     7517        assert (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound);
     7518      } else {
     7519        assert (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound);
     7520      }
     7521    }
     7522  }
     7523  int numberInBlock = block->numberInBlock_;
     7524  for (; i < numberInBlock; i++) {
     7525    int iColumn = column[i];
     7526    if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
     7527      assert (model->getColumnStatus(iColumn) == ClpSimplex::basic ||
     7528              model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
     7529    assert (lookup[iColumn] == i);
     7530  }
    62597531#endif
    62607532}
     
    62647536void
    62657537ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
    6266                                  const double * pi,
    6267                                  CoinIndexedVector * output) const
    6268 {
    6269      int numberNonZero = 0;
    6270      int * index = output->getIndices();
    6271      double * array = output->denseVector();
    6272      double zeroTolerance = model->zeroTolerance();
    6273      double value = 0.0;
    6274      CoinBigIndex j;
    6275      int numberOdd = block_->startIndices_;
    6276      if (numberOdd) {
    6277           // A) as probably long may be worth unrolling
    6278           CoinBigIndex end = start_[1];
    6279           for (j = start_[0]; j < end; j++) {
    6280                int iRow = row_[j];
    6281                value += pi[iRow] * element_[j];
    6282           }
    6283           int iColumn;
    6284           // int jColumn=column_[0];
    6285 
    6286           for (iColumn = 0; iColumn < numberOdd - 1; iColumn++) {
    6287                CoinBigIndex start = end;
    6288                end = start_[iColumn+2];
    6289                if (fabs(value) > zeroTolerance) {
    6290                     array[numberNonZero] = value;
    6291                     index[numberNonZero++] = column_[iColumn];
    6292                     //index[numberNonZero++]=jColumn;
    6293                }
    6294                // jColumn = column_[iColumn+1];
    6295                value = 0.0;
    6296                //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
    6297                for (j = start; j < end; j++) {
    6298                     int iRow = row_[j];
    6299                     value += pi[iRow] * element_[j];
    6300                }
    6301                //}
    6302           }
    6303           if (fabs(value) > zeroTolerance) {
    6304                array[numberNonZero] = value;
    6305                index[numberNonZero++] = column_[iColumn];
    6306                //index[numberNonZero++]=jColumn;
    6307           }
    6308      }
    6309      for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
    6310           // B) Can sort so just do nonbasic (and nonfixed)
    6311           // C) Can do two at a time (if so put odd one into start_)
    6312           // D) can use switch
    6313           blockStruct * block = block_ + iBlock;
    6314           //int numberPrice = block->numberInBlock_;
    6315           int numberPrice = block->numberPrice_;
    6316           int nel = block->numberElements_;
    6317           int * row = row_ + block->startElements_;
    6318           double * element = element_ + block->startElements_;
    6319           int * column = column_ + block->startIndices_;
    6320 #if 0
    6321           // two at a time
    6322           if ((numberPrice & 1) != 0) {
    6323                double value = 0.0;
    6324                int nel2 = nel;
    6325                for (; nel2; nel2--) {
    6326                     int iRow = *row++;
    6327                     value += pi[iRow] * (*element++);
    6328                }
    6329                if (fabs(value) > zeroTolerance) {
    6330                     array[numberNonZero] = value;
    6331                     index[numberNonZero++] = *column;
    6332                }
    6333                column++;
    6334           }
    6335           numberPrice = numberPrice >> 1;
    6336           switch ((nel % 2)) {
    6337                // 2 k +0
    6338           case 0:
    6339                for (; numberPrice; numberPrice--) {
    6340                     double value0 = 0.0;
    6341                     double value1 = 0.0;
    6342                     int nel2 = nel;
    6343                     for (; nel2; nel2--) {
    6344                          int iRow0 = *row;
    6345                          int iRow1 = *(row + nel);
    6346                          row++;
    6347                          double element0 = *element;
    6348                          double element1 = *(element + nel);
    6349                          element++;
    6350                          value0 += pi[iRow0] * element0;
    6351                          value1 += pi[iRow1] * element1;
    6352                     }
    6353                     row += nel;
    6354                     element += nel;
    6355                     if (fabs(value0) > zeroTolerance) {
    6356                          array[numberNonZero] = value0;
    6357                          index[numberNonZero++] = *column;
    6358                     }
    6359                     column++;
    6360                     if (fabs(value1) > zeroTolerance) {
    6361                          array[numberNonZero] = value1;
    6362                          index[numberNonZero++] = *column;
    6363                     }
    6364                     column++;
    6365                }
    6366                break;
    6367                // 2 k +1
    6368           case 1:
    6369                for (; numberPrice; numberPrice--) {
    6370                     double value0;
    6371                     double value1;
    6372                     int nel2 = nel - 1;
    6373                     {
    6374                          int iRow0 = row[0];
    6375                          int iRow1 = row[nel];
    6376                          double pi0 = pi[iRow0];
    6377                          double pi1 = pi[iRow1];
    6378                          value0 = pi0 * element[0];
    6379                          value1 = pi1 * element[nel];
    6380                          row++;
    6381                          element++;
    6382                     }
    6383                     for (; nel2; nel2--) {
    6384                          int iRow0 = *row;
    6385                          int iRow1 = *(row + nel);
    6386                          row++;
    6387                          double element0 = *element;
    6388                          double element1 = *(element + nel);
    6389                          element++;
    6390                          value0 += pi[iRow0] * element0;
    6391                          value1 += pi[iRow1] * element1;
    6392                     }
    6393                     row += nel;
    6394                     element += nel;
    6395                     if (fabs(value0) > zeroTolerance) {
    6396                          array[numberNonZero] = value0;
    6397                          index[numberNonZero++] = *column;
    6398                     }
    6399                     column++;
    6400                     if (fabs(value1) > zeroTolerance) {
    6401                          array[numberNonZero] = value1;
    6402                          index[numberNonZero++] = *column;
    6403                     }
    6404                     column++;
    6405                }
    6406                break;
    6407           }
     7538                                 const double * pi,
     7539                                 CoinIndexedVector * output) const
     7540{
     7541  int numberNonZero = 0;
     7542  int * index = output->getIndices();
     7543  double * array = output->denseVector();
     7544  double zeroTolerance = model->zeroTolerance();
     7545  double value = 0.0;
     7546  CoinBigIndex j;
     7547  int numberOdd = block_->startIndices_;
     7548  if (numberOdd) {
     7549    // A) as probably long may be worth unrolling
     7550    CoinBigIndex end = start_[1];
     7551    for (j = start_[0]; j < end; j++) {
     7552      int iRow = row_[j];
     7553      value += pi[iRow] * element_[j];
     7554    }
     7555    int iColumn;
     7556    // int jColumn=column_[0];
     7557   
     7558    for (iColumn = 0; iColumn < numberOdd - 1; iColumn++) {
     7559      CoinBigIndex start = end;
     7560      end = start_[iColumn+2];
     7561      if (fabs(value) > zeroTolerance) {
     7562        array[numberNonZero] = value;
     7563        index[numberNonZero++] = column_[iColumn];
     7564        //index[numberNonZero++]=jColumn;
     7565      }
     7566      // jColumn = column_[iColumn+1];
     7567      value = 0.0;
     7568      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
     7569      for (j = start; j < end; j++) {
     7570        int iRow = row_[j];
     7571        value += pi[iRow] * element_[j];
     7572      }
     7573      //}
     7574    }
     7575    if (fabs(value) > zeroTolerance) {
     7576      array[numberNonZero] = value;
     7577      index[numberNonZero++] = column_[iColumn];
     7578      //index[numberNonZero++]=jColumn;
     7579    }
     7580  }
     7581  for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
     7582    // C) Can do two at a time (if so put odd one into start_)
     7583    // D) can use switch
     7584    blockStruct * block = block_ + iBlock;
     7585    //int numberPrice = block->numberInBlock_;
     7586    int numberPrice = block->firstBasic_;
     7587    int nel = block->numberElements_;
     7588    int * row = row_ + block->startElements_;
     7589    double * element = element_ + block->startElements_;
     7590    int * column = column_ + block->startIndices_;
     7591    int nBlock=numberPrice>>COIN_AVX2_SHIFT;
     7592    numberPrice -= roundDown(numberPrice);
     7593    for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7594      for (int j=0; j < COIN_AVX2; j++) {
     7595        double value = 0.0;
     7596        for (int i=0; i < nel; i++) {
     7597          int iRow = row[i*COIN_AVX2];
     7598          value += pi[iRow] * element[i*COIN_AVX2];
     7599        }
     7600#if COIN_AVX2>1
     7601        row++;
     7602        element++;
    64087603#else
    6409           for (; numberPrice; numberPrice--) {
    6410                double value = 0.0;
    6411                int nel2 = nel;
    6412                for (; nel2; nel2--) {
    6413                     int iRow = *row++;
    6414                     value += pi[iRow] * (*element++);
    6415                }
    6416                if (fabs(value) > zeroTolerance) {
    6417                     array[numberNonZero] = value;
    6418                     index[numberNonZero++] = *column;
    6419                }
    6420                column++;
    6421           }
    6422 #endif
    6423      }
    6424      output->setNumElements(numberNonZero);
     7604        row+=nel;
     7605        element+=nel;
     7606#endif
     7607        if (fabs(value) > zeroTolerance) {
     7608          array[numberNonZero] = value;
     7609          index[numberNonZero++] = *column;
     7610        }
     7611        column++;
     7612      }
     7613#if COIN_AVX2>1
     7614      row += (nel-1)*COIN_AVX2;
     7615      element += (nel-1)*COIN_AVX2;
     7616      assert (row == row_ + block->startElements_ + nel*COIN_AVX2*(jBlock+1));
     7617#endif
     7618    }
     7619    // last lot
     7620    for (int j=0; j < numberPrice; j++) {
     7621      double value = 0.0;
     7622      for (int i=0; i < nel; i++) {
     7623        int iRow = row[i*COIN_AVX2];
     7624        value += pi[iRow] * element[i*COIN_AVX2];
     7625      }
     7626#if COIN_AVX2>1
     7627      row++;
     7628      element++;
     7629#else
     7630      row+=nel;
     7631      element+=nel;
     7632#endif
     7633      if (fabs(value) > zeroTolerance) {
     7634        array[numberNonZero] = value;
     7635        index[numberNonZero++] = *column;
     7636      }
     7637      column++;
     7638    }
     7639  }
     7640  output->setNumElements(numberNonZero);
     7641}
     7642static void
     7643transposeTimes3Bit2Odd(clpTempInfo & info)
     7644{
     7645  double zeroTolerance = info.tolerance;
     7646  double dualTolerance = -info.dualTolerance; // sign swapped
     7647  double * COIN_RESTRICT reducedCost = info.reducedCost;
     7648  double * COIN_RESTRICT weights = info.solution;
     7649  const unsigned int * reference =
     7650    reinterpret_cast<const unsigned int *>(info.upper);
     7651  const unsigned char * status = info.status;
     7652  const CoinBigIndex * starts = info.start;
     7653  const int * row = info.row;
     7654  const int * column = info.which;
     7655  const double * element = info.element;
     7656  double scaleFactor = info.theta;
     7657  const double * pi = info.cost;
     7658  const double * piWeight = info.lower;
     7659  double referenceIn = info.upperTheta;
     7660  double devex = info.changeObj;
     7661  assert (scaleFactor);
     7662  int bestSequence = info.numberAdded;
     7663  double bestRatio = info.bestPossible;
     7664#define NO_CHANGE_MULTIPLIER 1
     7665#define INVERSE_MULTIPLIER 1.0/NO_CHANGE_MULTIPLIER
     7666#if NO_CHANGE_MULTIPLIER != 1
     7667  double bestRatio2 = NO_CHANGE_MULTIPLIER*bestRatio;
     7668#else
     7669#define bestRatio2 bestRatio
     7670#endif
     7671  CoinBigIndex end = starts[0];
     7672  int numberOdd = info.numberToDo;
     7673  for (int jColumn = 0; jColumn < numberOdd; jColumn++) {
     7674    CoinBigIndex start = end;
     7675    CoinBigIndex j;
     7676    int iColumn = column[jColumn];
     7677    end = starts[jColumn+1];
     7678    double value = 0.0;
     7679    if ((status[iColumn]&7) != 1) {
     7680      for (j = start; j < end; j++) {
     7681        int iRow = row[j];
     7682        value -= pi[iRow] * element[j];
     7683      }
     7684      if (fabs(value) > zeroTolerance) {
     7685        // and do other array
     7686        double modification = 0.0;
     7687        for (j = start; j < end; j++) {
     7688          int iRow = row[j];
     7689          modification += piWeight[iRow] * element[j];
     7690        }
     7691        double thisWeight = weights[iColumn];
     7692        double pivot = value * scaleFactor;
     7693        double pivotSquared = pivot * pivot;
     7694        thisWeight += pivotSquared * devex + pivot * modification;
     7695        if (thisWeight < DEVEX_TRY_NORM) {
     7696          if (referenceIn < 0.0) {
     7697            // steepest
     7698            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     7699          } else {
     7700            // exact
     7701            thisWeight = referenceIn * pivotSquared;
     7702            if (reference(iColumn))
     7703              thisWeight += 1.0;
     7704            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     7705          }
     7706        }
     7707        weights[iColumn] = thisWeight;
     7708        value = reducedCost[iColumn]-value;
     7709        reducedCost[iColumn] = value;
     7710        unsigned char thisStatus=status[iColumn]&7;
     7711        if (thisStatus==3) {
     7712        } else if ((thisStatus&1)!=0) {
     7713          // basic or fixed
     7714          value=0.0;
     7715        } else if (thisStatus==2) {
     7716          value=-value;
     7717        } else {
     7718          // free or superbasic
     7719          if (fabs(value) > -FREE_ACCEPT * dualTolerance) {
     7720            // we are going to bias towards free (but only if reasonable)
     7721            value = -fabs(value)*FREE_BIAS;
     7722          } else {
     7723            value=0.0;
     7724          }
     7725        }
     7726        if (value < dualTolerance) {
     7727          value *= value;
     7728          if (value>bestRatio*weights[iColumn]) {
     7729            bestSequence = iColumn;
     7730            bestRatio = value/weights[iColumn];
     7731#if NO_CHANGE_MULTIPLIER != 1
     7732            bestRatio2 = bestRatio * NO_CHANGE_MULTIPLIER;
     7733#endif
     7734          }
     7735        }
     7736      } else {
     7737        value = reducedCost[iColumn];
     7738        unsigned char thisStatus=status[iColumn]&7;
     7739        if (thisStatus==3) {
     7740        } else if ((thisStatus&1)!=0) {
     7741          // basic or fixed
     7742          value=0.0;
     7743        } else if (thisStatus==2) {
     7744          value=-value;
     7745        } else {
     7746          // free or superbasic
     7747          if (fabs(value) > -FREE_ACCEPT * dualTolerance) {
     7748            // we are going to bias towards free (but only if reasonable)
     7749            value = -fabs(value)*FREE_BIAS;
     7750          } else {
     7751            value=0.0;
     7752          }
     7753        }
     7754        if (value < dualTolerance) {
     7755          value *= value;
     7756          if (value>bestRatio2*weights[iColumn]) {
     7757            bestSequence = iColumn;
     7758            bestRatio2 = value/weights[iColumn];
     7759#if NO_CHANGE_MULTIPLIER != 1
     7760            bestRatio = bestRatio2 * INVERSE_MULTIPLIER;
     7761#endif
     7762          }
     7763        }
     7764      }
     7765    }
     7766  }
     7767  info.numberAdded = bestSequence;
     7768  info.bestPossible = bestRatio;
     7769}
     7770#include <immintrin.h>
     7771static void
     7772transposeTimes3Bit2(clpTempInfo & info)
     7773{
     7774  double zeroTolerance = info.tolerance;
     7775  double dualTolerance = -info.dualTolerance; // sign swapped
     7776  double * COIN_RESTRICT reducedCost = info.reducedCost;
     7777  double * COIN_RESTRICT weights = info.solution;
     7778  double * COIN_RESTRICT work = info.work;
     7779  double * COIN_RESTRICT work2 = work+COIN_AVX2_CHUNK;
     7780  const unsigned int * reference =
     7781    reinterpret_cast<const unsigned int *>(info.upper);
     7782  const blockStruct * blocks =
     7783    reinterpret_cast<const blockStruct *>(info.pivotVariable);
     7784  const unsigned char * status = info.status;
     7785  const int * rowBlock = info.row;
     7786  const int * columnBlock = info.which;
     7787  const double * elementBlock = info.element;
     7788  double scaleFactor = info.theta;
     7789  const double * pi = info.cost;
     7790  const double * piWeight = info.lower;
     7791  double referenceIn = info.upperTheta;
     7792  double devex = info.changeObj;
     7793  assert (scaleFactor);
     7794  int bestSequence = info.numberAdded;
     7795  double bestRatio = info.bestPossible;
     7796#if NO_CHANGE_MULTIPLIER != 1
     7797  double bestRatio2 = NO_CHANGE_MULTIPLIER*bestRatio;
     7798#endif
     7799#define INCLUDE_MATRIX3_PRICING 1
     7800  int firstBlock=info.startColumn;
     7801  int lastBlock=info.numberToDo;
     7802  //double * COIN_RESTRICT tempArray =
     7803  //const_cast<double *>(elementBlock+info.numberColumns);
     7804  //double * COIN_RESTRICT tempArray2 = tempArray + roundUp(maxBlockSize);
     7805  for (int iBlock = firstBlock; iBlock < lastBlock; iBlock++) {
     7806    const blockStruct * block = blocks + iBlock;
     7807    int numberPrice = block->firstBasic_;
     7808    int nel = block->numberElements_;
     7809    const int * row = rowBlock + block->startElements_;
     7810    const double * element = elementBlock + block->startElements_;
     7811    const int * column = columnBlock + block->startIndices_;
     7812#if COIN_AVX2 > 1
     7813    int numberDo = roundDown(numberPrice);
     7814#ifdef HASWELL
     7815    assert(COIN_AVX2==4);
     7816    __m256d zero = _mm256_setzero_pd();
     7817    for (int kColumn=0;kColumn<numberDo;kColumn+=COIN_AVX2_CHUNK) {
     7818      int endBlock=CoinMin(COIN_AVX2_CHUNK,numberDo-kColumn);
     7819      for(int i=0;i<endBlock;i+=COIN_AVX2) {
     7820        __m256d arrayX = _mm256_setzero_pd();
     7821        __m256d tempX = _mm256_setzero_pd();
     7822        for (int j=0;j<nel;j++) {
     7823          __m128i rows = _mm_loadu_si128((const __m128i *)row);
     7824          __m256d elements = _mm256_load_pd(element);
     7825          __m256d pis = _mm256_i32gather_pd(pi,rows,8);
     7826          __m256d piWeights = _mm256_i32gather_pd(piWeight,rows,8);
     7827          // should be better way using sub
     7828          arrayX = _mm256_fmadd_pd(pis,elements,arrayX);
     7829          tempX = _mm256_fmadd_pd(piWeights,elements,tempX);
     7830          row+=4;
     7831          element+=4;
     7832        }
     7833        arrayX = _mm256_sub_pd(zero,arrayX);
     7834        _mm256_store_pd(work+i,tempX);
     7835        _mm256_store_pd(work2+i,arrayX);
     7836      }
     7837      for (int i=0;i<endBlock;i++) {
     7838        double value=work2[i];
     7839        double modification = work[i];
     7840        // common coding
     7841#include "ClpPackedMatrix.hpp"
     7842      }
     7843    }
     7844#else
     7845    for (int kColumn=0;kColumn<numberDo;kColumn+=COIN_AVX2_CHUNK) {
     7846      int n=0;
     7847      int nBlock=CoinMin(COIN_AVX2_CHUNK,numberPrice-kColumn)>>COIN_AVX2_SHIFT;
     7848      //int nBlock = numberPrice>>COIN_AVX2_SHIFT;
     7849      for (int jBlock=0;jBlock<nBlock;jBlock++) {
     7850        for (int j=0; j < COIN_AVX2; j++) {
     7851          double value = 0.0;
     7852          double modification = 0.0;
     7853          for (int i=0; i < nel; i++) {
     7854            //assert (element+i*COIN_AVX2<tempArray);
     7855            int iRow = row[i*COIN_AVX2];
     7856            value -= pi[iRow] * element[i*COIN_AVX2];
     7857            modification += piWeight[iRow] * element[i*COIN_AVX2];
     7858          }
     7859          work[n]=modification;
     7860          work2[n++] = value;
     7861          row++;
     7862          element++;
     7863        }
     7864        row += (nel-1)*COIN_AVX2;
     7865        element += (nel-1)*COIN_AVX2;
     7866        //assert (row == rowBlock + block->startElements_ + nel*COIN_AVX2*(jBlock+1));
     7867      }
     7868      for (int j=0;j<n;j++) {
     7869        double value=work2[j];
     7870        double modification = work[j];
     7871        // common coding
     7872#include "ClpPackedMatrix.hpp"
     7873      }
     7874    }
     7875#endif
     7876    // last lot
     7877    for (int j=numberDo; j < numberPrice; j++) {
     7878      double value = 0.0;
     7879      double modification = 0.0;
     7880      for (int i=0; i < nel; i++) {
     7881        int iRow = row[i*COIN_AVX2];
     7882        value -= pi[iRow] * element[i*COIN_AVX2];
     7883        modification += piWeight[iRow] * element[i*COIN_AVX2];
     7884      }
     7885      // common coding
     7886      #include "ClpPackedMatrix.hpp"
     7887      row++;
     7888      element++;
     7889    }
     7890#else // COIN_AVX2 == 1
     7891    for (int j=0; j < numberPrice; j++) {
     7892      double value = 0.0;
     7893      double modification = 0.0;
     7894      for (int i=0; i < nel; i++) {
     7895        //assert (element+i*COIN_AVX2<tempArray);
     7896        int iRow = row[i*COIN_AVX2];
     7897        value -= pi[iRow] * element[i*COIN_AVX2];
     7898        modification += piWeight[iRow] * element[i*COIN_AVX2];
     7899      }
     7900      // common coding
     7901      #include "ClpPackedMatrix.hpp"
     7902      row+=nel;
     7903      element+=nel;
     7904    }
     7905    int n = numberPrice;
     7906#endif
     7907  }
     7908  info.numberAdded = bestSequence;
     7909  info.bestPossible = bestRatio;
     7910}
     7911static void
     7912transposeTimes3BitSlacks(clpTempInfo & info)
     7913{
     7914  double dualTolerance = info.dualTolerance;
     7915  double * COIN_RESTRICT reducedCost = info.reducedCost;
     7916  double * COIN_RESTRICT weights = info.solution;
     7917  const int * columnBlock = info.which;
     7918  const blockStruct * blocks =
     7919    reinterpret_cast<const blockStruct *>(info.pivotVariable);
     7920  //const unsigned char * status = info.status;
     7921  int bestSequence = info.numberAdded;
     7922  double bestRatio = info.bestPossible;
     7923  int iBlock=info.startColumn;
     7924  assert (info.numberToDo==iBlock+1);
     7925  const blockStruct * block = blocks + iBlock;
     7926  int first = 0;
     7927  int last = block->firstAtUpper_;
     7928  const int * column = columnBlock + block->startIndices_;
     7929  for (int j=first;j<last;j++) {
     7930    int iColumn=*column;
     7931    column++;
     7932    double value = reducedCost[iColumn];
     7933    // free or superbasic
     7934    if (fabs(value) > FREE_ACCEPT * dualTolerance) {
     7935      // we are going to bias towards free (but only if reasonable)
     7936      value = -fabs(value)*FREE_BIAS;
     7937      value *= value;
     7938      if (value>bestRatio*weights[iColumn]) {
     7939        bestSequence = iColumn;
     7940        bestRatio = value/weights[iColumn];
     7941      }
     7942    }
     7943  }
     7944  first = last;
     7945  last = block->firstAtLower_;
     7946  for (int j=first;j<last;j++) {
     7947    int iColumn=*column;
     7948    column++;
     7949    double value = reducedCost[iColumn];
     7950    // at upper
     7951    if (value > dualTolerance) {
     7952      value *= value;
     7953      if (value>bestRatio*weights[iColumn]) {
     7954        bestSequence = iColumn;
     7955        bestRatio = value/weights[iColumn];
     7956      }
     7957    }
     7958  }
     7959  first = last;
     7960  last = block->firstBasic_;
     7961  dualTolerance = - dualTolerance; //flip sign
     7962  for (int j=first;j<last;j++) {
     7963    int iColumn=*column;
     7964    column++;
     7965    double value = reducedCost[iColumn];
     7966    // at lower
     7967    if (value < dualTolerance) {
     7968      value *= value;
     7969      if (value>bestRatio*weights[iColumn]) {
     7970        bestSequence = iColumn;
     7971        bestRatio = value/weights[iColumn];
     7972      }
     7973    }
     7974  }
     7975  info.numberAdded = bestSequence;
     7976  info.bestPossible = bestRatio;
    64257977}
    64267978// Updates two arrays for steepest
    64277979void
    64287980ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
    6429                                   const double * pi, CoinIndexedVector * output,
    6430                                   const double * piWeight,
    6431                                   double referenceIn, double devex,
    6432                                   // Array for exact devex to say what is in reference framework
    6433                                   unsigned int * reference,
    6434                                   double * weights, double scaleFactor)
    6435 {
    6436      int numberNonZero = 0;
    6437      int * index = output->getIndices();
    6438      double * array = output->denseVector();
    6439      double zeroTolerance = model->zeroTolerance();
    6440      double value = 0.0;
    6441      bool killDjs = (scaleFactor == 0.0);
    6442      if (!scaleFactor)
    6443           scaleFactor = 1.0;
    6444      int numberOdd = block_->startIndices_;
    6445      int iColumn;
    6446      CoinBigIndex end = start_[0];
    6447      for (iColumn = 0; iColumn < numberOdd; iColumn++) {
    6448           CoinBigIndex start = end;
    6449           CoinBigIndex j;
    6450           int jColumn = column_[iColumn];
    6451           end = start_[iColumn+1];
    6452           value = 0.0;
    6453           if (model->getColumnStatus(jColumn) != ClpSimplex::basic) {
    6454                for (j = start; j < end; j++) {
    6455                     int iRow = row_[j];
    6456                     value -= pi[iRow] * element_[j];
    6457                }
    6458                if (fabs(value) > zeroTolerance) {
    6459                     // and do other array
    6460                     double modification = 0.0;
    6461                     for (j = start; j < end; j++) {
    6462                          int iRow = row_[j];
    6463                          modification += piWeight[iRow] * element_[j];
    6464                     }
    6465                     double thisWeight = weights[jColumn];
    6466                     double pivot = value * scaleFactor;
    6467                     double pivotSquared = pivot * pivot;
    6468                     thisWeight += pivotSquared * devex + pivot * modification;
    6469                     if (thisWeight < DEVEX_TRY_NORM) {
    6470                          if (referenceIn < 0.0) {
    6471                               // steepest
    6472                               thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    6473                          } else {
    6474                               // exact
    6475                               thisWeight = referenceIn * pivotSquared;
    6476                               if (reference(jColumn))
    6477                                    thisWeight += 1.0;
    6478                               thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    6479                          }
    6480                     }
    6481                     weights[jColumn] = thisWeight;
    6482                     if (!killDjs) {
    6483                          array[numberNonZero] = value;
    6484                          index[numberNonZero++] = jColumn;
    6485                     }
    6486                }
    6487           }
    6488      }
    6489      for (int iBlock = 0; iBlock < numberBlocks_; iBlock++) {
    6490           // B) Can sort so just do nonbasic (and nonfixed)
    6491           // C) Can do two at a time (if so put odd one into start_)
    6492           // D) can use switch
    6493           blockStruct * block = block_ + iBlock;
    6494           //int numberPrice = block->numberInBlock_;
    6495           int numberPrice = block->numberPrice_;
    6496           int nel = block->numberElements_;
    6497           int * row = row_ + block->startElements_;
    6498           double * element = element_ + block->startElements_;
    6499           int * column = column_ + block->startIndices_;
    6500           for (; numberPrice; numberPrice--) {
    6501                double value = 0.0;
    6502                int nel2 = nel;
    6503                for (; nel2; nel2--) {
    6504                     int iRow = *row++;
    6505                     value -= pi[iRow] * (*element++);
    6506                }
    6507                if (fabs(value) > zeroTolerance) {
    6508                     int jColumn = *column;
    6509                     // back to beginning
    6510                     row -= nel;
    6511                     element -= nel;
    6512                     // and do other array
    6513                     double modification = 0.0;
    6514                     nel2 = nel;
    6515                     for (; nel2; nel2--) {
    6516                          int iRow = *row++;
    6517                          modification += piWeight[iRow] * (*element++);
    6518                     }
    6519                     double thisWeight = weights[jColumn];
    6520                     double pivot = value * scaleFactor;
    6521                     double pivotSquared = pivot * pivot;
    6522                     thisWeight += pivotSquared * devex + pivot * modification;
    6523                     if (thisWeight < DEVEX_TRY_NORM) {
    6524                          if (referenceIn < 0.0) {
    6525                               // steepest
    6526                               thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    6527                          } else {
    6528                               // exact
    6529                               thisWeight = referenceIn * pivotSquared;
    6530                               if (reference(jColumn))
    6531                                    thisWeight += 1.0;
    6532                               thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    6533                          }
    6534                     }
    6535                     weights[jColumn] = thisWeight;
    6536                     if (!killDjs) {
    6537                          array[numberNonZero] = value;
    6538                          index[numberNonZero++] = jColumn;
    6539                     }
    6540                }
    6541                column++;
    6542           }
    6543      }
    6544      output->setNumElements(numberNonZero);
    6545      output->setPackedMode(true);
     7981                                  const double * COIN_RESTRICT pi,
     7982                                  CoinIndexedVector * output,
     7983                                  const double * COIN_RESTRICT piWeight,
     7984                                  double * COIN_RESTRICT infeas,
     7985                                  double * COIN_RESTRICT reducedCost,
     7986                                  double referenceIn, double devex,
     7987                                  // Array for exact devex to say what is in reference framework
     7988                                  unsigned int * COIN_RESTRICT reference,
     7989                                  double * COIN_RESTRICT weights, double scaleFactor)
     7990{
     7991  /*
     7992    This is designed for meaty problems so probably should also go to exact devex
     7993    also could do sequenceIn_ choice here
     7994    could do ones at lb separately - ?? masking - ? masking for weights
     7995   */
     7996  double zeroTolerance = model->zeroTolerance();
     7997  double dualTolerance = model->currentDualTolerance();
     7998  // we can't really trust infeasibilities if there is dual error
     7999  // this coding has to mimic coding in checkDualSolution
     8000  double error = CoinMin(1.0e-2, model->largestDualError());
     8001  // allow tolerance at least slightly bigger than standard
     8002  dualTolerance = dualTolerance +  error;
     8003  int numberOdd = block_->startIndices_;
     8004  assert (scaleFactor);
     8005#if ABOCA_LITE
     8006  //int nChunks=numberChunks_;
     8007#define ODD_INFO 2*ABOCA_LITE
     8008#else
     8009#define ODD_INFO 1
     8010  //int nChunks=1;
     8011#endif
     8012  clpTempInfo info[ODD_INFO+2]={};
     8013  unsigned char * status = model->statusArray();
     8014  // fill in info
     8015  for (int i=0;i<ODD_INFO+2;i++) {
     8016    //memset(info+i,0,sizeof(clpTempInfo));
     8017    info[i].tolerance=zeroTolerance;
     8018    info[i].dualTolerance=dualTolerance;
     8019    info[i].reducedCost=reducedCost;
     8020    info[i].infeas=infeas;
     8021    info[i].solution=weights;
     8022    info[i].upper=reinterpret_cast<const double *>(reference);
     8023    info[i].pivotVariable=reinterpret_cast<const int *>(block_);
     8024    info[i].status=status;
     8025    info[i].start=start_;
     8026    info[i].row=row_;
     8027    info[i].which=column_;
     8028    info[i].element=element_;
     8029    info[i].theta=scaleFactor;
     8030    info[i].cost=pi;
     8031    info[i].lower=piWeight;
     8032    info[i].upperTheta=referenceIn;
     8033    info[i].changeObj=devex;
     8034    info[i].numberAdded=-1;
     8035  }
     8036  //assert (ODD_INFO==1);
     8037  info[ODD_INFO].startColumn=0;
     8038  info[ODD_INFO].numberToDo=numberOdd;
     8039  //
     8040  double * temporary = temporary_->array();
     8041#if ABOCA_LITE
     8042  for (int iBlock=0;iBlock<numberChunks_;iBlock++) {
     8043    info[iBlock].startColumn=endChunk_[iBlock];
     8044    info[iBlock].numberToDo=endChunk_[iBlock+1];
     8045    info[iBlock].work=temporary+2*COIN_AVX2_CHUNK*iBlock;
     8046  }
     8047#else
     8048  info[0].startColumn=0;
     8049  info[0].numberToDo=numberBlocks_;
     8050  info[0].work=temporary;
     8051#endif
     8052  info[ODD_INFO+1].startColumn=numberBlocks_;
     8053  info[ODD_INFO+1].numberToDo=numberBlocks_+1;
     8054#if ABOCA_LITE
     8055  if (abcState()>1) {
     8056    cilk_spawn transposeTimes3Bit2Odd(info[ODD_INFO]);
     8057    for (int iBlock=0;iBlock<numberChunks_;iBlock++) {
     8058      cilk_spawn transposeTimes3Bit2(info[iBlock]);
     8059    }
     8060    transposeTimes3BitSlacks(info[ODD_INFO+1]);
     8061    cilk_sync;
     8062  } else {
     8063#endif
     8064    transposeTimes3Bit2Odd(info[ODD_INFO]);
     8065#if ABOCA_LITE
     8066    for (int iBlock=0;iBlock<numberChunks_;iBlock++) {
     8067      transposeTimes3Bit2(info[iBlock]);
     8068    }
     8069#else
     8070    transposeTimes3Bit2(info[0]);
     8071#endif
     8072    transposeTimes3BitSlacks(info[ODD_INFO+1]);
     8073#if ABOCA_LITE
     8074  }
     8075#endif 
     8076  int sequenceIn=-1;
     8077  double best=0.0;
     8078  for (int iBlock=0;iBlock<ODD_INFO+2;iBlock++) {
     8079    if (info[iBlock].bestPossible>best) {
     8080      best=info[iBlock].bestPossible;
     8081      sequenceIn=info[iBlock].numberAdded;
     8082    }
     8083  }
     8084  // Make sure sequenceOut not eligible
     8085  int sequenceOut = model->sequenceOut();
     8086  double saveOutDj;
     8087  if (sequenceOut>=0) {
     8088    saveOutDj=reducedCost[sequenceOut];
     8089    unsigned char thisStatus=status[sequenceOut]&7;
     8090    assert (thisStatus!=0&&thisStatus!=4);
     8091    if (thisStatus!=2)
     8092      reducedCost[sequenceOut]=COIN_DBL_MAX;
     8093    else
     8094      reducedCost[sequenceOut]=-COIN_DBL_MAX;
     8095  }
     8096  if (sequenceIn>=0) {
     8097    // check not flagged
     8098    if (model->flagged(sequenceIn)||sequenceIn==sequenceOut) {
     8099      int number = model->numberRows() + model->numberColumns();
     8100      const unsigned char * COIN_RESTRICT status =
     8101        model->statusArray();
     8102      sequenceIn = -2;
     8103      double bestRatio=0.0;
     8104      for (int iSequence=0;iSequence<number;iSequence++) {
     8105        unsigned char thisStatus=status[iSequence]&7;
     8106        double value = reducedCost[iSequence];
     8107        if (thisStatus==3) {
     8108        } else if ((thisStatus&1)!=0) {
     8109          // basic or fixed
     8110          value=0.0;
     8111        } else if (thisStatus==2) {
     8112          value=-value;
     8113        } else {
     8114          // free or superbasic
     8115          if (fabs(value) > FREE_ACCEPT * -dualTolerance) {
     8116            // we are going to bias towards free (but only if reasonable)
     8117            value = -fabs(value)*FREE_BIAS;
     8118          } else {
     8119            value=0.0;
     8120          }
     8121        }
     8122        if (value < dualTolerance) {
     8123          value *= value;
     8124          if (value>bestRatio*weights[iSequence]) {
     8125            if (!model->flagged(iSequence)) {
     8126              sequenceIn = iSequence;
     8127              bestRatio = value/weights[iSequence];
     8128            }
     8129          }
     8130        }
     8131      }
     8132    }
     8133  }
     8134  if (sequenceOut>=0)
     8135    reducedCost[sequenceOut]=saveOutDj;
     8136  model->spareIntArray_[3]=sequenceIn;
     8137}
     8138 /*
     8139  type - 1 redo infeasible, 2 choose sequenceIn, 3 both
     8140  returns sequenceIn (or -1) for type 2
     8141 */
     8142int
     8143ClpPackedMatrix3::redoInfeasibilities(const ClpSimplex * model,
     8144                                      ClpPrimalColumnSteepest * pivotChoose,
     8145                                      int type)
     8146{
     8147  CoinIndexedVector * infeasible = pivotChoose->infeasible();
     8148  // change to use blocks
     8149  double tolerance = model->currentDualTolerance();
     8150  // we can't really trust infeasibilities if there is dual error
     8151  // this coding has to mimic coding in checkDualSolution
     8152  double error = CoinMin(1.0e-2, model->largestDualError());
     8153  // allow tolerance at least slightly bigger than standard
     8154  tolerance = tolerance +  error;
     8155  // reverse sign so test is cleaner
     8156  tolerance = - tolerance;
     8157  int number = model->numberRows() + model->numberColumns();
     8158
     8159  const double * COIN_RESTRICT reducedCost = model->djRegion();
     8160  const unsigned char * COIN_RESTRICT status =
     8161    model->statusArray();
     8162  const double * COIN_RESTRICT weights = pivotChoose->weights();
     8163  int sequenceIn=-1;
     8164  double bestRatio=0.0;
     8165  switch (type) {
     8166  case 1:
     8167    infeasible->clear();
     8168    for (int iSequence=0;iSequence<number;iSequence++) {
     8169      unsigned char thisStatus=status[iSequence]&7;
     8170      double value = reducedCost[iSequence];
     8171      if (thisStatus==3) {
     8172      } else if ((thisStatus&1)!=0) {
     8173        // basic or fixed
     8174        value=0.0;
     8175      } else if (thisStatus==2) {
     8176        value=-value;
     8177      } else {
     8178        // free or superbasic
     8179        if (fabs(value) > FREE_ACCEPT * -tolerance) {
     8180          // we are going to bias towards free (but only if reasonable)
     8181          value = -fabs(value)*FREE_BIAS;
     8182        } else {
     8183          value=0.0;
     8184        }
     8185      }
     8186      if (value < tolerance) {
     8187        // store square in list
     8188        infeasible->quickAdd(iSequence, value * value);
     8189      }
     8190    }
     8191    break;
     8192  case 2:
     8193    infeasible->clear(); // temp for debug
     8194    for (int iSequence=0;iSequence<number;iSequence++) {
     8195      unsigned char thisStatus=status[iSequence]&7;
     8196      double value = reducedCost[iSequence];
     8197      if (thisStatus==3) {
     8198      } else if ((thisStatus&1)!=0) {
     8199        // basic or fixed
     8200        value=0.0;
     8201      } else if (thisStatus==2) {
     8202        value=-value;
     8203      } else {
     8204        // free or superbasic
     8205        if (fabs(value) > FREE_ACCEPT * -tolerance) {
     8206          // we are going to bias towards free (but only if reasonable)
     8207          value = -fabs(value)*FREE_BIAS;
     8208        } else {
     8209          value=0.0;
     8210        }
     8211      }
     8212      if (value < tolerance) {
     8213        value *= value;
     8214        if (value>bestRatio*weights[iSequence]) {
     8215          sequenceIn = iSequence;
     8216          bestRatio = value/weights[iSequence];
     8217        }
     8218      }
     8219    }
     8220    break;
     8221  case 3:
     8222    infeasible->clear();
     8223    for (int iSequence=0;iSequence<number;iSequence++) {
     8224      unsigned char thisStatus=status[iSequence]&7;
     8225      double value = reducedCost[iSequence];
     8226      if (thisStatus==3) {
     8227      } else if ((thisStatus&1)!=0) {
     8228        // basic or fixed
     8229        value=0.0;
     8230      } else if (thisStatus==2) {
     8231        value=-value;
     8232      } else {
     8233        // free or superbasic
     8234        if (fabs(value) > FREE_ACCEPT * -tolerance) {
     8235          // we are going to bias towards free (but only if reasonable)
     8236          value = -fabs(value)*FREE_BIAS;
     8237        } else {
     8238          value=0.0;
     8239        }
     8240      }
     8241      if (value < tolerance) {
     8242        value *= value;
     8243        // store square in list
     8244        infeasible->quickAdd(iSequence, value);
     8245        if (value>bestRatio*weights[iSequence]) {
     8246          sequenceIn = iSequence;
     8247          bestRatio = value/weights[iSequence];
     8248        }
     8249      }
     8250    }
     8251    break;
     8252  }
     8253  if (sequenceIn>=0) {
     8254    // check not flagged
     8255    if (model->flagged(sequenceIn)) {
     8256      sequenceIn = -1;
     8257      for (int iSequence=0;iSequence<number;iSequence++) {
     8258        unsigned char thisStatus=status[iSequence]&7;
     8259        double value = reducedCost[iSequence];
     8260        if (thisStatus==3) {
     8261        } else if ((thisStatus&1)!=0) {
     8262          // basic or fixed
     8263          value=0.0;
     8264        } else if (thisStatus==2) {
     8265          value=-value;
     8266        } else {
     8267          // free or superbasic
     8268          if (fabs(value) > FREE_ACCEPT * -tolerance) {
     8269            // we are going to bias towards free (but only if reasonable)
     8270            value = -fabs(value)*FREE_BIAS;
     8271          } else {
     8272            value=0.0;
     8273          }
     8274        }
     8275        if (value < tolerance) {
     8276          value *= value;
     8277          if (value>bestRatio*weights[iSequence]) {
     8278            if (!model->flagged(iSequence)) {
     8279              sequenceIn = iSequence;
     8280              bestRatio = value/weights[iSequence];
     8281            }
     8282          }
     8283        }
     8284      }
     8285    }
     8286  }
     8287  return sequenceIn;
    65468288}
    65478289#if COIN_LONG_WORK
    65488290// For long double versions
    65498291void
    6550 ClpPackedMatrix::times(CoinWorkDouble scalar,
    6551                       const CoinWorkDouble * x, CoinWorkDouble * y) const
    6552 {
    6553      int iRow, iColumn;
    6554      // get matrix data pointers
    6555      const int * row = matrix_->getIndices();
    6556      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    6557      const double * elementByColumn = matrix_->getElements();
    6558      //memset(y,0,matrix_->getNumRows()*sizeof(double));
    6559      assert (((flags_ & 2) != 0) == matrix_->hasGaps());
    6560      if (!(flags_ & 2)) {
    6561           for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    6562                CoinBigIndex j;
    6563                CoinWorkDouble value = x[iColumn];
    6564                if (value) {
    6565                     CoinBigIndex start = columnStart[iColumn];
    6566                     CoinBigIndex end = columnStart[iColumn+1];
    6567                     value *= scalar;
    6568                     for (j = start; j < end; j++) {
    6569                          iRow = row[j];
    6570                          y[iRow] += value * elementByColumn[j];
    6571                     }
    6572                }
    6573           }
    6574      } else {
    6575           const int * columnLength = matrix_->getVectorLengths();
    6576           for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    6577                CoinBigIndex j;
    6578                CoinWorkDouble value = x[iColumn];
    6579                if (value) {
    6580                     CoinBigIndex start = columnStart[iColumn];
    6581                     CoinBigIndex end = start + columnLength[iColumn];
    6582                     value *= scalar;
    6583                     for (j = start; j < end; j++) {
    6584                          iRow = row[j];
    6585                          y[iRow] += value * elementByColumn[j];
    6586                     }
    6587                }
    6588           }
    6589      }
     8292  ClpPackedMatrix::times(CoinWorkDouble scalar,
     8293                        const CoinWorkDouble * x, CoinWorkDouble * y) const
     8294{
     8295  int iRow, iColumn;
     8296  // get matrix data pointers
     8297  const int * row = matrix_->getIndices();
     8298  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     8299  const double * elementByColumn = matrix_->getElements();
     8300  //memset(y,0,matrix_->getNumRows()*sizeof(double));
     8301  assert (((flags_ & 2) != 0) == matrix_->hasGaps());
     8302  if (!(flags_ & 2)) {
     8303    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     8304      CoinBigIndex j;
     8305      CoinWorkDouble value = x[iColumn];
     8306      if (value) {
     8307        CoinBigIndex start = columnStart[iColumn];
     8308        CoinBigIndex end = columnStart[iColumn+1];
     8309        value *= scalar;
     8310        for (j = start; j < end; j++) {
     8311          iRow = row[j];
     8312          y[iRow] += value * elementByColumn[j];
     8313        }
     8314      }
     8315    }
     8316  } else {
     8317    const int * columnLength = matrix_->getVectorLengths();
     8318    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     8319      CoinBigIndex j;
     8320      CoinWorkDouble value = x[iColumn];
     8321      if (value) {
     8322        CoinBigIndex start = columnStart[iColumn];
     8323        CoinBigIndex end = start + columnLength[iColumn];
     8324        value *= scalar;
     8325        for (j = start; j < end; j++) {
     8326          iRow = row[j];
     8327          y[iRow] += value * elementByColumn[j];
     8328        }
     8329      }
     8330    }
     8331  }
    65908332}
    65918333void
    65928334ClpPackedMatrix::transposeTimes(CoinWorkDouble scalar,
    6593                                 const CoinWorkDouble * x, CoinWorkDouble * y) const
    6594 {
    6595      int iColumn;
    6596      // get matrix data pointers
    6597      const int * row = matrix_->getIndices();
    6598      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    6599      const double * elementByColumn = matrix_->getElements();
    6600      if (!(flags_ & 2)) {
    6601           if (scalar == -1.0) {
    6602                CoinBigIndex start = columnStart[0];
    6603                for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    6604                     CoinBigIndex j;
    6605                     CoinBigIndex next = columnStart[iColumn+1];
    6606                     CoinWorkDouble value = y[iColumn];
    6607                     for (j = start; j < next; j++) {
    6608                          int jRow = row[j];
    6609                          value -= x[jRow] * elementByColumn[j];
    6610                     }
    6611                     start = next;
    6612                     y[iColumn] = value;
    6613                }
    6614           } else {
    6615                CoinBigIndex start = columnStart[0];
    6616                for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    6617                     CoinBigIndex j;
    6618                     CoinBigIndex next = columnStart[iColumn+1];
    6619                     CoinWorkDouble value = 0.0;
    6620                     for (j = start; j < next; j++) {
    6621                          int jRow = row[j];
    6622                          value += x[jRow] * elementByColumn[j];
    6623                     }
    6624                     start = next;
    6625                     y[iColumn] += value * scalar;
    6626                }
    6627           }
    6628      } else {
    6629           const int * columnLength = matrix_->getVectorLengths();
    6630           for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
    6631                CoinBigIndex j;
    6632                CoinWorkDouble value = 0.0;
    6633                CoinBigIndex start = columnStart[iColumn];
    6634                CoinBigIndex end = start + columnLength[iColumn];
    6635                for (j = start; j < end; j++) {
    6636                     int jRow = row[j];
    6637                     value += x[jRow] * elementByColumn[j];
    6638                }
    6639                y[iColumn] += value * scalar;
    6640           }
    6641      }
     8335                                const CoinWorkDouble * x, CoinWorkDouble * y) const
     8336{
     8337  int iColumn;
     8338  // get matrix data pointers
     8339  const int * row = matrix_->getIndices();
     8340  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     8341  const double * elementByColumn = matrix_->getElements();
     8342  if (!(flags_ & 2)) {
     8343    if (scalar == -1.0) {
     8344      CoinBigIndex start = columnStart[0];
     8345      for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     8346        CoinBigIndex j;
     8347        CoinBigIndex next = columnStart[iColumn+1];
     8348        CoinWorkDouble value = y[iColumn];
     8349        for (j = start; j < next; j++) {
     8350          int jRow = row[j];
     8351          value -= x[jRow] * elementByColumn[j];
     8352        }
     8353        start = next;
     8354        y[iColumn] = value;
     8355      }
     8356    } else {
     8357      CoinBigIndex start = columnStart[0];
     8358      for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     8359        CoinBigIndex j;
     8360        CoinBigIndex next = columnStart[iColumn+1];
     8361        CoinWorkDouble value = 0.0;
     8362        for (j = start; j < next; j++) {
     8363          int jRow = row[j];
     8364          value += x[jRow] * elementByColumn[j];
     8365        }
     8366        start = next;
     8367        y[iColumn] += value * scalar;
     8368      }
     8369    }
     8370  } else {
     8371    const int * columnLength = matrix_->getVectorLengths();
     8372    for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) {
     8373      CoinBigIndex j;
     8374      CoinWorkDouble value = 0.0;
     8375      CoinBigIndex start = columnStart[iColumn];
     8376      CoinBigIndex end = start + columnLength[iColumn];
     8377      for (j = start; j < end; j++) {
     8378        int jRow = row[j];
     8379        value += x[jRow] * elementByColumn[j];
     8380      }
     8381      y[iColumn] += value * scalar;
     8382    }
     8383  }
    66428384}
    66438385#endif
     
    66458387#undef reference
    66468388#endif
     8389 
     8390 
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r2078 r2235  
    1010
    1111#include "ClpMatrixBase.hpp"
     12#include "ClpPrimalColumnSteepest.hpp"
    1213
    1314/** This implements CoinPackedMatrix as derived from ClpMatrixBase.
     
    1920class ClpPackedMatrix2;
    2021class ClpPackedMatrix3;
     22class CoinDoubleArrayWithLength;
    2123class ClpPackedMatrix : public ClpMatrixBase {
    2224
     
    123125     /** Creates scales for column copy (rowCopy in model may be modified)
    124126         returns non-zero if no scaling done */
    125      virtual int scale(ClpModel * model, const ClpSimplex * baseModel = NULL) const ;
     127     virtual int scale(ClpModel * model, ClpSimplex * simplex = NULL) const ;
    126128     /** Scales rowCopy if column copy scaled
    127129         Only called if scales already exist */
     
    263265     virtual bool canCombine(const ClpSimplex * model,
    264266                             const CoinIndexedVector * pi) const;
    265      /// Updates two arrays for steepest
    266      virtual void transposeTimes2(const ClpSimplex * model,
    267                                   const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    268                                   const CoinIndexedVector * pi2,
    269                                   CoinIndexedVector * spare,
     267     /** Updates two arrays for steepest and does devex weights
     268         Returns nonzero if updates reduced cost and infeas -
     269         new infeas in dj1 */
     270     virtual int transposeTimes2(const ClpSimplex * model,
     271                                 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
     272                                 const CoinIndexedVector * pi2,
     273                                 CoinIndexedVector * spare,
     274                                 double * infeas, double * reducedCost,
    270275                                  double referenceIn, double devex,
    271276                                  // Array for exact devex to say what is in reference framework
     
    561566};
    562567typedef struct {
    563      CoinBigIndex startElements_; // point to data
    564      int startIndices_; // point to column_
    565      int numberInBlock_;
    566      int numberPrice_; // at beginning
    567      int numberElements_; // number elements per column
     568  CoinBigIndex startElements_; // point to data
     569  CoinBigIndex startRows_; // point to data later
     570  int startIndices_; // point to column_
     571  int numberInBlock_;
     572  int numberScan_; // i.e. miss out basic and fixed
     573  /* order is -
     574     free or superbasic
     575     at upper
     576     at lower
     577     fixed or basic */
     578  int firstAtUpper_;
     579  int firstAtLower_;
     580  int firstBasic_; // or fixed
     581  int numberElements_; // number elements per column
     582  int numberOnes_; // later
    568583} blockStruct;
    569584class ClpPackedMatrix3 {
     
    582597                          const double * pi, CoinIndexedVector * dj1,
    583598                          const double * piWeight,
     599                          double * COIN_RESTRICT infeas,
     600                          double * COIN_RESTRICT reducedCost,
    584601                          double referenceIn, double devex,
    585602                          // Array for exact devex to say what is in reference framework
     
    612629     void swapOne(const ClpSimplex * model, const ClpPackedMatrix * matrix,
    613630                  int iColumn);
     631     /// Part of above
     632     void swapOne(int iBlock, int kA, int kB);
     633     /**
     634        type - 1 redo infeasible, 2 choose sequenceIn, 3 both
     635        returns sequenceIn (or -1) for type 2
     636     */
     637     int redoInfeasibilities(const ClpSimplex * model,
     638                             ClpPrimalColumnSteepest * pivotChoose,
     639                             int type);
     640  /// Get temporary array (aligned)
    614641     //@}
    615642
     
    623650     /// Number of columns
    624651     int numberColumns_;
     652     /// Number of columns including gaps
     653     int numberColumnsWithGaps_;
     654#if ABOCA_LITE
     655     /// Number of chunks
     656     int numberChunks_;
     657#endif
     658     /// Number of elements (including gaps)
     659     CoinBigIndex numberElements_;
     660     /// Maximum size of any block
     661     int maxBlockSize_;
    625662     /// Column indices and reverse lookup (within block)
    626663     int * column_;
    627      /// Starts for odd/long vectors
     664     /// Starts for odd/long vectors??
    628665     CoinBigIndex * start_;
    629666     /// Rows
     
    631668     /// Elements
    632669     double * element_;
     670     /// Temporary work area (aligned)
     671     CoinDoubleArrayWithLength * temporary_;
     672#if ABOCA_LITE
     673     /// Chunk ends (could have more than cpus)
     674     int endChunk_[2*ABOCA_LITE+1];
     675#endif
    633676     /// Blocks (ordinary start at 0 and go to first block)
    634677     blockStruct * block_;
    635678     //@}
    636679};
    637 
    638 #endif
     680#elif INCLUDE_MATRIX3_PRICING
     681      int iColumn=*column;
     682      column++;
     683      if (fabs(value) > zeroTolerance) {
     684        double thisWeight = weights[iColumn];
     685        double pivot = value * scaleFactor;
     686        double pivotSquared = pivot * pivot;
     687        thisWeight += pivotSquared * devex + pivot * modification;
     688        if (thisWeight < DEVEX_TRY_NORM) {
     689          if (referenceIn < 0.0) {
     690            // steepest
     691            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
     692          } else {
     693            // exact
     694            thisWeight = referenceIn * pivotSquared;
     695            if (reference(iColumn))
     696              thisWeight += 1.0;
     697            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
     698          }
     699        }
     700        // out basic or fixed
     701        weights[iColumn] = thisWeight;
     702        value = reducedCost[iColumn]-value;
     703        reducedCost[iColumn] = value;
     704        unsigned char thisStatus=status[iColumn]&7;
     705        assert (thisStatus!=0&&thisStatus!=4);
     706        if (thisStatus==3) {
     707          //} else if ((thisStatus&1)!=0) {
     708          // basic or fixed
     709          //value=0.0;
     710        } else {
     711          assert (thisStatus==2);
     712          value=-value;
     713        }
     714        if (value < dualTolerance) {
     715          value *= value;
     716          if (value>bestRatio*weights[iColumn]) {
     717            bestSequence = iColumn;
     718            bestRatio = value/weights[iColumn];
     719#if NO_CHANGE_MULTIPLIER != 1
     720            bestRatio2 = bestRatio * NO_CHANGE_MULTIPLIER;
     721#endif
     722          }
     723        }
     724      } else {
     725        // interesting - was faster without this?!
     726        value = reducedCost[iColumn];
     727        unsigned char thisStatus=status[iColumn]&7;
     728        assert (thisStatus!=0&&thisStatus!=4);
     729        if (thisStatus==3) {
     730        } else if ((thisStatus&1)!=0) {
     731          // basic or fixed
     732          value=0.0;
     733        } else {
     734          value=-value;
     735        }
     736        if (value < dualTolerance) {
     737          value *= value;
     738          if (value>bestRatio2*weights[iColumn]) {
     739            bestSequence = iColumn;
     740            bestRatio2 = value/weights[iColumn];
     741#if NO_CHANGE_MULTIPLIER != 1
     742            bestRatio = bestRatio2 * INVERSE_MULTIPLIER;
     743#endif
     744          }
     745        }
     746      }
     747#endif
  • trunk/Clp/src/ClpPlusMinusOneMatrix.cpp

    r2078 r2235  
    19261926// These have to match ClpPrimalColumnSteepest version
    19271927#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
    1928 // Updates two arrays for steepest
    1929 void
     1928/* Updates two arrays for steepest and does devex weights
     1929   Returns nonzero if updates reduced cost and infeas -
     1930   new infeas in dj1  - This does not at present*/
     1931int
    19301932ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex * model,
    19311933                                       const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    19321934                                       const CoinIndexedVector * pi2,
    19331935                                       CoinIndexedVector * spare,
     1936                                       double * , double * ,
    19341937                                       double referenceIn, double devex,
    19351938                                       // Array for exact devex to say what is in reference framework
     
    20712074     if (packed)
    20722075          dj1->setPackedMode(true);
     2076     return 0;
    20732077}
    20742078// Updates second array for steepest and does devex weights
  • trunk/Clp/src/ClpPlusMinusOneMatrix.hpp

    r2078 r2235  
    178178     virtual bool canCombine(const ClpSimplex * model,
    179179                             const CoinIndexedVector * pi) const;
    180      /// Updates two arrays for steepest
    181      virtual void transposeTimes2(const ClpSimplex * model,
    182                                   const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    183                                   const CoinIndexedVector * pi2,
    184                                   CoinIndexedVector * spare,
     180     /** Updates two arrays for steepest and does devex weights
     181         Returns nonzero if updates reduced cost and infeas -
     182         new infeas in dj1 */
     183     virtual int transposeTimes2(const ClpSimplex * model,
     184                                 const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
     185                                 const CoinIndexedVector * pi2,
     186                                 CoinIndexedVector * spare,
     187                                 double * infeas, double * reducedCost,
    185188                                  double referenceIn, double devex,
    186189                                  // Array for exact devex to say what is in reference framework
  • trunk/Clp/src/ClpPresolve.cpp

    r2186 r2235  
    10261026          }
    10271027          if (ifree) {
    1028             int fill_level=10;
     1028            int fill_level=CoinMax(10,prob->maxSubstLevel_);
    10291029            const CoinPresolveAction * lastAction = NULL;
    10301030            int iPass=4;
     
    13471347                         // do costed if Clp (at end as ruins rest of presolve)
    13481348                      possibleBreak;
     1349#ifndef CLP_MOVE_COSTS
    13491350                         paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
     1351#else
     1352                         double * fakeRowObjective=new double[prob->nrows_];
     1353                         memset(fakeRowObjective,0,prob->nrows_*sizeof(double));
     1354                         paction_ = slack_singleton_action::presolve(prob, paction_, fakeRowObjective);
     1355                         delete [] fakeRowObjective;
     1356#endif
    13501357                         stopLoop = true;
    13511358                    }
  • trunk/Clp/src/ClpPrimalColumnSteepest.cpp

    r2074 r2235  
    55
    66#include "CoinPragma.hpp"
     7#if ABOCA_LITE
     8// 1 is not owner of abcState_
     9#define ABCSTATE_LITE 1
     10#endif
     11//#define FAKE_CILK
    712
    813#include "ClpSimplex.hpp"
     
    1217#include "ClpNonLinearCost.hpp"
    1318#include "ClpMessage.hpp"
     19#include "ClpEventHandler.hpp"
     20#include "ClpPackedMatrix.hpp"
    1421#include "CoinHelperFunctions.hpp"
    1522#include <stdio.h>
     23#if ABOCA_LITE
     24#undef ALT_UPDATE_WEIGHTS
     25#define ALT_UPDATE_WEIGHTS 2
     26#endif
     27#if ALT_UPDATE_WEIGHTS==1
     28extern CoinIndexedVector * altVector[3];
     29#endif
     30static void debug1(int iSequence,double value,double weight)
     31{
     32  printf("xx %d inf %.20g wt %.20g\n",
     33         iSequence,value,weight);
     34}
    1635//#define CLP_DEBUG
    1736//#############################################################################
     
    3251       state_(-1),
    3352       mode_(mode),
     53       infeasibilitiesState_(0),
    3454       persistence_(normal),
    3555       numberSwitched_(0),
     
    4969     state_ = rhs.state_;
    5070     mode_ = rhs.mode_;
     71     infeasibilitiesState_ = rhs.infeasibilitiesState_;
    5172     persistence_ = rhs.persistence_;
    5273     numberSwitched_ = rhs.numberSwitched_;
     
    115136          state_ = rhs.state_;
    116137          mode_ = rhs.mode_;
     138          infeasibilitiesState_ = rhs.infeasibilitiesState_;
    117139          persistence_ = rhs.persistence_;
    118140          numberSwitched_ = rhs.numberSwitched_;
     
    160182#define TRY_NORM 1.0e-4
    161183#define ADD_ONE 1.0
     184static void
     185pivotColumnBit(clpTempInfo & info)
     186{
     187  double * COIN_RESTRICT weights = const_cast<double *>(info.lower);
     188  const unsigned char * COIN_RESTRICT status = info.status;
     189  double tolerance=info.tolerance;
     190  double bestDj=info.primalRatio;
     191  int bestSequence=-1;
     192  double * COIN_RESTRICT infeas = const_cast<double *>(info.infeas);
     193  const int * COIN_RESTRICT start = info.which;
     194  const int * COIN_RESTRICT index = info.index;
     195  for (int i = start[0]; i < start[1]; i++) {
     196    int iSequence = index[i];
     197    double value = infeas[iSequence];
     198    double weight = weights[iSequence];
     199    if (value > tolerance) {
     200      if (value > bestDj * weight) {
     201        // check flagged variable and correct dj
     202        if ((status[iSequence]&64)==0) {
     203          bestDj = value / weight;
     204          bestSequence = iSequence;
     205        }
     206      }
     207    }
     208  }
     209  info.primalRatio=bestDj;
     210  info.numberAdded=bestSequence;
     211}
    162212// Returns pivot column, -1 if none
    163213/*      The Packed CoinIndexedVector updates has cost updates - for normal LP
     
    287337          }
    288338     }
     339     int bestSequence = -1;
     340     model_->spareIntArray_[3] = -3;
    289341     if (switchType == 5) {
    290342          if (anyUpdates > 0) {
     
    297349               djsAndSteepest(updates, spareRow2,
    298350                              spareColumn1, spareColumn2);
     351               if (model_->spareIntArray_[3]>-2) {
     352                 bestSequence = model_->spareIntArray_[3];
     353                 infeasibilitiesState_=2;
     354               } else if (model_->spareIntArray_[3]==-2) {
     355                 redoInfeasibilities();
     356               }
    299357          } else {
    300358               // devex etc when can use dj
     
    323381          }
    324382     }
     383     // everything may have been done if vector copy
     384     if (infeasibilitiesState_ == 2) {
     385       // all done
     386       infeasibilitiesState_ = 1;
     387       model_->clpMatrix()->setSavedBestSequence(bestSequence);
     388       if (bestSequence >= 0)
     389         model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
     390       assert (sequenceOut != bestSequence);
     391       return bestSequence;
     392     } else if (infeasibilitiesState_ == 1) {
     393       // need to redo
     394       //infeasibilitiesState_ = 0;
     395       redoInfeasibilities();
     396     }
     397       
    325398#ifdef CLP_DEBUG
    326399     alternateWeights_->checkClear();
     
    440513     if (switchType < 4) {
    441514          if (switchType < 2 ) {
    442                numberWanted = number + 1;
     515               numberWanted = COIN_INT_MAX-1;
    443516          } else if (switchType == 2) {
    444517               numberWanted = CoinMax(2000, number / 8);
     
    469542
    470543     double bestDj = 1.0e-30;
    471      int bestSequence = -1;
     544     bestSequence = -1;
     545#ifdef CLP_USER_DRIVEN
     546     // could go parallel?
     547     model_->eventHandler()->event(ClpEventHandler::beforeChooseIncoming);
     548#endif
    472549
    473550     int i, iSequence;
     
    516593
    517594     int iPass;
    518      // Setup two passes
    519      int start[4];
    520      start[1] = number;
    521      start[2] = 0;
    522      double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
    523      start[0] = static_cast<int> (dstart);
    524      start[3] = start[0];
     595     // Setup two passes (unless all)
     596     if (mode_>1) {
     597       int start[4];
     598       start[1] = number;
     599       start[2] = 0;
     600       double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
     601       start[0] = static_cast<int> (dstart);
     602       start[3] = start[0];
     603       for (iPass = 0; iPass < 2; iPass++) {
     604         int end = start[2*iPass+1];
     605         if (switchType < 5) {
     606           for (i = start[2*iPass]; i < end; i++) {
     607             iSequence = index[i];
     608             double value = infeas[iSequence];
     609             double weight = weights_[iSequence];
     610             if (value > tolerance) {
     611               //weight=1.0;
     612               if (value > bestDj * weight) {
     613                 // check flagged variable and correct dj
     614                 if (!model_->flagged(iSequence)) {
     615                   bestDj = value / weight;
     616                   bestSequence = iSequence;
     617                 } else {
     618                   // just to make sure we don't exit before got something
     619                   numberWanted++;
     620                 }
     621               }
     622               numberWanted--;
     623             }
     624             if (!numberWanted)
     625               break;
     626           }
     627         } else {
     628           // Dantzig
     629           for (i = start[2*iPass]; i < end; i++) {
     630             iSequence = index[i];
     631             double value = infeas[iSequence];
     632             if (value > tolerance) {
     633               if (value > bestDj) {
     634                 // check flagged variable and correct dj
     635                 if (!model_->flagged(iSequence)) {
     636                   bestDj = value;
     637                   bestSequence = iSequence;
     638                 } else {
     639                   // just to make sure we don't exit before got something
     640                   numberWanted++;
     641                 }
     642               }
     643               numberWanted--;
     644             }
     645             if (!numberWanted)
     646               break;
     647           }
     648         }
     649         if (!numberWanted)
     650           break;
     651       }
     652     } else {
     653#if ABOCA_LITE
     654       int numberThreads=CoinMax(abcState(),1);
     655#define ABOCA_LITE_MAX ABOCA_LITE
     656#else
     657       const int numberThreads=1;
     658#define ABOCA_LITE_MAX 1
     659#endif
     660       if (0) {
     661         int iSequence;
     662         double value;
     663         double weight;
     664         iSequence=34841;
     665         value = infeas[iSequence];
     666         weight = weights_[iSequence];
     667         debug1(iSequence,value,weight);
     668         iSequence=34845;
     669         value = infeas[iSequence];
     670         weight = weights_[iSequence];
     671         debug1(iSequence,value,weight);
     672       }
     673       clpTempInfo info[ABOCA_LITE_MAX];
     674       int start_lite[2*ABOCA_LITE_MAX];
     675       int chunk0 = (number+numberThreads-1)/numberThreads;
     676       int n0=0;
     677       for (int i=0;i<numberThreads;i++) {
     678         int * startX=start_lite+2*i;
     679         info[i].primalRatio = bestDj;
     680         info[i].lower = weights_;
     681         info[i].infeas = infeas;
     682         info[i].index = index;
     683         info[i].status=const_cast<unsigned char *>(model_->statusArray());
     684         info[i].which=startX;
     685         info[i].tolerance=tolerance;
     686         startX[0]=n0;
     687         startX[1]=CoinMin(n0+chunk0,number);
     688         n0 += chunk0;
     689       }
     690       if (numberThreads==1) {
     691         pivotColumnBit(info[0]);
     692       } else {
     693         for (int i=0;i<numberThreads;i++) {
     694           cilk_spawn pivotColumnBit(info[i]);
     695         }
     696         cilk_sync;
     697       }
     698       for (int i=0;i<numberThreads;i++) {
     699         double bestDjX = info[i].primalRatio;
     700         if (bestDjX>bestDj) {
     701           bestDj=bestDjX;
     702           bestSequence=info[i].numberAdded;
     703         }
     704       }
     705     }
    525706     //double largestWeight=0.0;
    526707     //double smallestWeight=1.0e100;
    527      for (iPass = 0; iPass < 2; iPass++) {
    528           int end = start[2*iPass+1];
    529           if (switchType < 5) {
    530                for (i = start[2*iPass]; i < end; i++) {
    531                     iSequence = index[i];
    532                     double value = infeas[iSequence];
    533                     double weight = weights_[iSequence];
    534                     if (value > tolerance) {
    535                          //weight=1.0;
    536                          if (value > bestDj * weight) {
    537                               // check flagged variable and correct dj
    538                               if (!model_->flagged(iSequence)) {
    539                                    bestDj = value / weight;
    540                                    bestSequence = iSequence;
    541                               } else {
    542                                    // just to make sure we don't exit before got something
    543                                    numberWanted++;
    544                               }
    545                          }
    546                          numberWanted--;
    547                     }
    548                     if (!numberWanted)
    549                          break;
    550                }
    551           } else {
    552                // Dantzig
    553                for (i = start[2*iPass]; i < end; i++) {
    554                     iSequence = index[i];
    555                     double value = infeas[iSequence];
    556                     if (value > tolerance) {
    557                          if (value > bestDj) {
    558                               // check flagged variable and correct dj
    559                               if (!model_->flagged(iSequence)) {
    560                                    bestDj = value;
    561                                    bestSequence = iSequence;
    562                               } else {
    563                                    // just to make sure we don't exit before got something
    564                                    numberWanted++;
    565                               }
    566                          }
    567                          numberWanted--;
    568                     }
    569                     if (!numberWanted)
    570                          break;
    571                }
    572           }
    573           if (!numberWanted)
    574                break;
    575      }
    576708     model_->clpMatrix()->setSavedBestSequence(bestSequence);
    577709     if (bestSequence >= 0)
     
    583715       printf("%d best %g\n",bestSequence,bestDj);*/
    584716
     717#ifdef CLP_USER_DRIVEN
     718     // swap when working
     719     struct {int bestSequence;double bestDj;} tempInfo;
     720     tempInfo.bestDj=bestDj;tempInfo.bestSequence=bestSequence;
     721     model_->eventHandler()->eventWithInfo(ClpEventHandler::afterChooseIncoming,
     722                                   reinterpret_cast<void *>(&tempInfo));
     723#endif
    585724#ifndef NDEBUG
    586725     if (bestSequence >= 0) {
     
    591730          }
    592731     }
     732#endif
     733#ifdef ALT_UPDATE_WEIGHTSz
     734     printf("weights");
     735     for (int i=0;i<model_->numberColumns()+model_->numberRows();i++) {
     736       if (model_->getColumnStatus(i)==ClpSimplex::isFixed||
     737           model_->getColumnStatus(i)==ClpSimplex::basic)
     738         weights_[i]=1.0;
     739       printf(" %g",weights_[i]);
     740       if ((i%10)==9)
     741         printf("\n");
     742     }
     743     printf("\n");
    593744#endif
    594745#if 0
     
    9691120     for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    9701121          if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    971                     !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1122                    model_->getStatus(iCheck) != ClpSimplex::isFixed)
    9721123               checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    9731124     }
     
    10021153     assert (updates->getIndices()[0] == pivotSequence_);
    10031154     pivotSequence_ = -1;
     1155#if 1
    10041156     //updates->scanAndPack();
    10051157     model_->factorization()->updateColumnTranspose(spareRow2, updates);
    10061158     //alternateWeights_->scanAndPack();
     1159#if ALT_UPDATE_WEIGHTS !=2
    10071160     model_->factorization()->updateColumnTranspose(spareRow2,
    10081161               alternateWeights_);
     1162#elif ALT_UPDATE_WEIGHTS==1
     1163     if (altVector[1]) {
     1164       int numberRows=model_->numberRows();
     1165       double * work1 = altVector[1]->denseVector();
     1166       double * worka = alternateWeights_->denseVector();
     1167       int iRow=-1;
     1168       double diff=1.0e-8;
     1169       for (int i=0;i<numberRows;i++) {
     1170         double dd=CoinMax(fabs(work1[i]),fabs(worka[i]));
     1171         double d=fabs(work1[i]-worka[i]);
     1172         if (dd>1.0e-6&&d>diff*dd) {
     1173           diff=d/dd;
     1174           iRow=i;
     1175         }
     1176       }
     1177       if (iRow>=0)
     1178         printf("largest difference of %g (%g,%g) on row %d\n",
     1179                diff,work1[iRow],worka[iRow],iRow);
     1180     }
     1181#endif
     1182#else
     1183     model_->factorization()->updateTwoColumnsTranspose(spareRow2, updates,alternateWeights_);
     1184#endif
    10091185     // and we can see if reference
    10101186     int sequenceIn = model_->sequenceIn();
     
    10281204     double * other = alternateWeights_->denseVector();
    10291205     int numberColumns = model_->numberColumns();
     1206     // if if (model_->clpMatrix()->canCombine(model_, pi1)) no need to index
    10301207     // rows
    10311208     reducedCost = model_->djRegion(0);
     
    10361213     updateBy = updates->denseVector();
    10371214     weight = weights_ + numberColumns;
     1215#ifdef CLP_USER_DRIVEN
     1216     model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
     1217#endif
    10381218
    10391219     for (j = 0; j < number; j++) {
     
    11691349          }
    11701350     }
     1351#ifdef CLP_USER_DRIVEN
     1352     // could go parallel?
     1353     //model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
     1354#endif
    11711355     // put row of tableau in rowArray and columnArray (packed)
    11721356     // get subset which have nonzero tableau elements
    1173      transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2,
    1174                      -scaleFactor);
     1357     int returnCode=transposeTimes2(updates, spareColumn1,
     1358                                    alternateWeights_, spareColumn2, spareRow2,
     1359                                    -scaleFactor);
    11751360     // zero updateBy
    11761361     CoinZeroN(updateBy, number);
     
    11781363     // columns
    11791364     assert (scaleFactor);
    1180      reducedCost = model_->djRegion(1);
    11811365     number = spareColumn1->getNumElements();
    11821366     index = spareColumn1->getIndices();
    11831367     updateBy = spareColumn1->denseVector();
    1184 
    1185      for (j = 0; j < number; j++) {
    1186           int iSequence = index[j];
    1187           double value = reducedCost[iSequence];
    1188           double value2 = updateBy[j];
    1189           updateBy[j] = 0.0;
    1190           value -= value2;
    1191           reducedCost[iSequence] = value;
    1192           ClpSimplex::Status status = model_->getStatus(iSequence);
    1193 
    1194           switch(status) {
    1195 
    1196           case ClpSimplex::basic:
    1197           case ClpSimplex::isFixed:
    1198                break;
    1199           case ClpSimplex::isFree:
    1200           case ClpSimplex::superBasic:
    1201                if (fabs(value) > FREE_ACCEPT * tolerance) {
    1202                     // we are going to bias towards free (but only if reasonable)
    1203                     value *= FREE_BIAS;
    1204                     // store square in list
    1205                     if (infeas[iSequence])
    1206                          infeas[iSequence] = value * value; // already there
    1207                     else
    1208                          infeasible_->quickAdd(iSequence, value * value);
    1209                } else {
    1210                     infeasible_->zero(iSequence);
    1211                }
    1212                break;
    1213           case ClpSimplex::atUpperBound:
    1214                if (value > tolerance) {
    1215                     // store square in list
    1216                     if (infeas[iSequence])
    1217                          infeas[iSequence] = value * value; // already there
    1218                     else
    1219                          infeasible_->quickAdd(iSequence, value * value);
    1220                } else {
    1221                     infeasible_->zero(iSequence);
    1222                }
    1223                break;
    1224           case ClpSimplex::atLowerBound:
    1225                if (value < -tolerance) {
    1226                     // store square in list
    1227                     if (infeas[iSequence])
    1228                          infeas[iSequence] = value * value; // already there
    1229                     else
    1230                          infeasible_->quickAdd(iSequence, value * value);
    1231                } else {
    1232                     infeasible_->zero(iSequence);
    1233                }
    1234           }
     1368     if (returnCode != 2 && infeasibilitiesState_) {
     1369       //spareColumn1->clear();
     1370       redoInfeasibilities();
     1371     }
     1372     if (returnCode == 1) {
     1373       // most work already done
     1374       for (j = 0; j < number; j++) {
     1375         int iSequence = index[j];
     1376         double value = updateBy[j];
     1377         if (value) {
     1378           updateBy[j]=0.0;
     1379           infeasible_->quickAdd(iSequence, value);
     1380         } else {
     1381           infeasible_->zero(iSequence);
     1382         }
     1383       }
     1384     } else if (returnCode==0) {
     1385       reducedCost = model_->djRegion(1);
     1386
     1387       for (j = 0; j < number; j++) {
     1388         int iSequence = index[j];
     1389         double value = reducedCost[iSequence];
     1390         double value2 = updateBy[j];
     1391         updateBy[j] = 0.0;
     1392         value -= value2;
     1393         reducedCost[iSequence] = value;
     1394         ClpSimplex::Status status = model_->getStatus(iSequence);
     1395         
     1396         switch(status) {
     1397           
     1398         case ClpSimplex::basic:
     1399         case ClpSimplex::isFixed:
     1400           break;
     1401         case ClpSimplex::isFree:
     1402         case ClpSimplex::superBasic:
     1403           if (fabs(value) > FREE_ACCEPT * tolerance) {
     1404             // we are going to bias towards free (but only if reasonable)
     1405             value *= FREE_BIAS;
     1406             // store square in list
     1407             if (infeas[iSequence])
     1408               infeas[iSequence] = value * value; // already there
     1409             else
     1410               infeasible_->quickAdd(iSequence, value * value);
     1411           } else {
     1412             infeasible_->zero(iSequence);
     1413           }
     1414           break;
     1415         case ClpSimplex::atUpperBound:
     1416           if (value > tolerance) {
     1417             // store square in list
     1418             if (infeas[iSequence])
     1419               infeas[iSequence] = value * value; // already there
     1420             else
     1421               infeasible_->quickAdd(iSequence, value * value);
     1422           } else {
     1423             infeasible_->zero(iSequence);
     1424           }
     1425           break;
     1426         case ClpSimplex::atLowerBound:
     1427           if (value < -tolerance) {
     1428             // store square in list
     1429             if (infeas[iSequence])
     1430               infeas[iSequence] = value * value; // already there
     1431             else
     1432               infeasible_->quickAdd(iSequence, value * value);
     1433           } else {
     1434             infeasible_->zero(iSequence);
     1435           }
     1436         }
     1437       }
     1438     } else {
     1439       //assert(!number);
    12351440     }
    12361441     // restore outgoing weight
    1237      if (sequenceOut >= 0)
    1238           weights_[sequenceOut] = outgoingWeight;
     1442     if (sequenceOut >= 0) {
     1443       //#define GROW_REFERENCE
     1444#ifdef GROW_REFERENCE
     1445       if (!reference(sequenceOut)) {
     1446         outgoingWeight += 1.0;
     1447         setReference(sequenceOut,true);
     1448       }
     1449#endif
     1450       weights_[sequenceOut] = outgoingWeight;
     1451       //if (model_->getStatus(sequenceOut) != ClpSimplex::basic &&
     1452       //  model_->getStatus(sequenceOut) != ClpSimplex::isFixed)
     1453       //checkAccuracy(sequenceOut, 1.0e-1, updates, spareRow2);
     1454     }
    12391455     // make sure infeasibility on incoming is 0.0
    12401456     infeasible_->zero(sequenceIn);
     
    12491465     for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    12501466          if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1251                     !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1467                    model_->getStatus(iCheck) != ClpSimplex::isFixed)
    12521468               checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    12531469     }
     
    14601676          for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    14611677               if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1462                          !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1678                         model_->getStatus(iCheck) != ClpSimplex::isFixed)
    14631679                    checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    14641680          }
     
    14961712     model_->clpMatrix()->transposeTimes(model_, -1.0,
    14971713                                         updates, spareColumn2, spareColumn1);
     1714#ifdef CLP_USER_DRIVEN
     1715     model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
     1716#endif
    14981717     // normal
    14991718     for (iSection = 0; iSection < 2; iSection++) {
     
    16391858          weight = weights_ + numberColumns;
    16401859          if (needSubset) {
    1641                // now update weight update array
    1642                model_->factorization()->updateColumnTranspose(spareRow2,
    1643                          alternateWeights_);
     1860#if ALT_UPDATE_WEIGHTS !=2
     1861            model_->factorization()->updateColumnTranspose(spareRow2,
     1862               alternateWeights_);
     1863#elif ALT_UPDATE_WEIGHTS==1
     1864     if (altVector[1]) {
     1865       int numberRows=model_->numberRows();
     1866       double * work1 = altVector[1]->denseVector();
     1867       double * worka = alternateWeights_->denseVector();
     1868       int iRow=-1;
     1869       double diff=1.0e-8;
     1870       for (int i=0;i<numberRows;i++) {
     1871         double dd=CoinMax(fabs(work1[i]),fabs(worka[i]));
     1872         double d=fabs(work1[i]-worka[i]);
     1873         if (dd>1.0e-6&&d>diff*dd) {
     1874           diff=d/dd;
     1875           iRow=i;
     1876         }
     1877       }
     1878       if (iRow>=0)
     1879         printf("largest2 difference of %g (%g,%g) on row %d\n",
     1880                diff,work1[iRow],worka[iRow],iRow);
     1881     }
     1882#endif
    16441883               // do alternateWeights_ here so can scale
    16451884               for (j = 0; j < number; j++) {
     
    17281967          for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    17291968               if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1730                          !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1969                         model_->getStatus(iCheck) != ClpSimplex::isFixed)
    17311970                    checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    17321971          }
     
    17371976}
    17381977// Updates two arrays for steepest
    1739 void
     1978int
    17401979ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    17411980          const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
     
    17541993          referenceIn = -1.0;
    17551994     }
     1995     int returnCode=0;
    17561996     if (model_->clpMatrix()->canCombine(model_, pi1)) {
    1757           // put row of tableau in rowArray and columnArray
    1758           model_->clpMatrix()->transposeTimes2(model_, pi1, dj1, pi2, spare, referenceIn, devex_,
     1997       double * infeas = scaleFactor ? infeasible_->denseVector() : NULL;
     1998       // put row of tableau in rowArray and columnArray
     1999       returnCode=model_->clpMatrix()->transposeTimes2(model_, pi1,
     2000                                                       dj1, pi2, spare,
     2001                                                       infeas,
     2002                                                       model_->djRegion(1),
     2003                                               referenceIn, devex_,
    17592004                                               reference_,
    17602005                                               weights_, scaleFactor);
     2006       if (model_->spareIntArray_[3]>-2)
     2007         returnCode=2;
    17612008     } else {
    17622009          // put row of tableau in rowArray and columnArray
     
    18112058     }
    18122059     dj2->setNumElements(0);
     2060     return returnCode;
    18132061}
    18142062// Update weights for Devex
     
    19122160     for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    19132161          if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1914                     !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     2162                    model_->getStatus(iCheck) != ClpSimplex::isFixed)
    19152163               checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    19162164     }
     
    19812229     // now update weight update array
    19822230     //alternateWeights_->scanAndPack();
     2231#if ALT_UPDATE_WEIGHTS !=2
    19832232     model_->factorization()->updateColumnTranspose(spareRow2,
    1984                alternateWeights_);
     2233                                                    alternateWeights_);
     2234#elif ALT_UPDATE_WEIGHTS==1
     2235     if (altVector[1]) {
     2236       int numberRows=model_->numberRows();
     2237       double * work1 = altVector[1]->denseVector();
     2238       double * worka = alternateWeights_->denseVector();
     2239       int iRow=-1;
     2240       double diff=1.0e-8;
     2241       for (int i=0;i<numberRows;i++) {
     2242         double dd=CoinMax(fabs(work1[i]),fabs(worka[i]));
     2243         double d=fabs(work1[i]-worka[i]);
     2244         if (dd>1.0e-6&&d>diff*dd) {
     2245           diff=d/dd;
     2246           iRow=i;
     2247         }
     2248       }
     2249       if (iRow>=0)
     2250         printf("largest3 difference of %g (%g,%g) on row %d\n",
     2251                diff,work1[iRow],worka[iRow],iRow);
     2252     }
     2253#endif
    19852254     // get subset which have nonzero tableau elements
    19862255     model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
     
    20582327     for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    20592328          if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    2060                     !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     2329                    model_->getStatus(iCheck) != ClpSimplex::isFixed)
    20612330               checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    20622331     }
     
    26202889               model_->factorization()->updateColumnTranspose(spareRow2,
    26212890                         alternateWeights_);
     2891#ifdef ALT_UPDATE_WEIGHTS
     2892               abort();
     2893#endif
    26222894               for (j = 0; j < number; j++) {
    26232895                    int iSequence = index[j];
     
    27272999          for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    27283000               if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    2729                          !model_->getStatus(iCheck) != ClpSimplex::isFixed)
     3001                         model_->getStatus(iCheck) != ClpSimplex::isFixed)
    27303002                    checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    27313003          }
     
    28613133     }
    28623134}
     3135void
     3136ClpPrimalColumnSteepest::redoInfeasibilities()
     3137{
     3138  double * COIN_RESTRICT infeas = infeasible_->denseVector();
     3139  int * COIN_RESTRICT index = infeasible_->getIndices();
     3140  double tolerance = model_->currentDualTolerance();
     3141  // we can't really trust infeasibilities if there is dual error
     3142  // this coding has to mimic coding in checkDualSolution
     3143  double error = CoinMin(1.0e-2, model_->largestDualError());
     3144  // allow tolerance at least slightly bigger than standard
     3145  tolerance = tolerance +  error;
     3146  // reverse sign so test is cleaner
     3147  tolerance = - tolerance;
     3148  int number = model_->numberRows() + model_->numberColumns();
     3149  int numberNonZero=0;
     3150  const double * COIN_RESTRICT reducedCost = model_->djRegion();
     3151  const unsigned char * COIN_RESTRICT status =
     3152    model_->statusArray();
     3153  for (int iSequence=0;iSequence<number;iSequence++) {
     3154    unsigned char thisStatus=status[iSequence]&7;
     3155    double value = reducedCost[iSequence];
     3156    infeas[iSequence]=0.0;
     3157    if (thisStatus==3) {
     3158    } else if ((thisStatus&1)!=0) {
     3159      // basic or fixed
     3160      value=0.0;
     3161    } else if (thisStatus==2) {
     3162      value=-value;
     3163    } else {
     3164      // free or superbasic
     3165      if (fabs(value) > FREE_ACCEPT * -tolerance) {
     3166        // we are going to bias towards free (but only if reasonable)
     3167        value = -fabs(value)*FREE_BIAS;
     3168      } else {
     3169        value=0.0;
     3170      }
     3171    }
     3172    if (value < tolerance) {
     3173      // store square in list
     3174      infeas[iSequence] = value * value;
     3175      index[numberNonZero++] = iSequence;
     3176    }
     3177  }
     3178  infeasible_->setNumElements(numberNonZero);
     3179  infeasibilitiesState_ = 0;
     3180}
    28633181/*
    28643182   1) before factorization
     
    28933211     bool doInfeasibilities = true;
    28943212     if (mode == 1) {
     3213#if ABOCA_LITE
     3214       int numberThreads=abcState();
     3215       if (numberThreads&&mode_>1)
     3216         mode_=0; // force exact
     3217#endif
    28953218          if(weights_) {
    28963219               // Check if size has changed
     
    29003223                    //alternateWeights_->clear();
    29013224                    if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
     3225#if ALT_UPDATE_WEIGHTS !=2
    29023226                         // save pivot order
    29033227                         CoinMemcpyN(pivotVariable,
    29043228                                     numberRows, alternateWeights_->getIndices());
     3229#endif
    29053230                         // change from pivot row number to sequence number
    29063231                         pivotSequence_ = pivotVariable[pivotSequence_];
     
    30023327          if (mode != 3) {
    30033328               if (pivotSequence_ >= 0) {
     3329#if ALT_UPDATE_WEIGHTS !=2
    30043330                    // restore pivot row
    30053331                    int iRow;
     
    30363362                         temp[iPivot] = 0.0;
    30373363                    }
     3364#else
     3365                    for (int iRow = 0; iRow < numberRows; iRow++) {
     3366                         int iPivot = pivotVariable[iRow];
     3367                         if (iPivot == pivotSequence_) {
     3368                           pivotSequence_ = iRow;
     3369                           break;
     3370                         }
     3371                    }
     3372#endif
    30383373               } else {
    30393374                    // Just clean up
     
    30473382          if(!doInfeasibilities)
    30483383               return; // don't disturb infeasibilities
    3049           infeasible_->clear();
     3384          double * COIN_RESTRICT infeas = infeasible_->denseVector();
     3385          int * COIN_RESTRICT index = infeasible_->getIndices();
     3386          int numberNonZero = 0;
     3387          infeasibilitiesState_ = 0;
    30503388          double tolerance = model_->currentDualTolerance();
    30513389          int number = model_->numberRows() + model_->numberColumns();
     
    30593397
    30603398          if (!model_->nonLinearCost()->lookBothWays()) {
     3399            const unsigned char * COIN_RESTRICT status =
     3400              model_->statusArray();
    30613401#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
    30623402               for (iSequence = 0; iSequence < number; iSequence++) {
    30633403                    double value = reducedCost[iSequence];
    3064                     ClpSimplex::Status status = model_->getStatus(iSequence);
    3065 
    3066                     switch(status) {
    3067 
    3068                     case ClpSimplex::basic:
    3069                     case ClpSimplex::isFixed:
    3070                          break;
    3071                     case ClpSimplex::isFree:
    3072                     case ClpSimplex::superBasic:
    3073                          if (fabs(value) > FREE_ACCEPT * tolerance) {
    3074                               // we are going to bias towards free (but only if reasonable)
    3075                               value *= FREE_BIAS;
    3076                               // store square in list
    3077                               infeasible_->quickAdd(iSequence, value * value);
    3078                          }
    3079                          break;
    3080                     case ClpSimplex::atUpperBound:
    3081                          if (value > tolerance) {
    3082                               infeasible_->quickAdd(iSequence, value * value);
    3083                          }
    3084                          break;
    3085                     case ClpSimplex::atLowerBound:
    3086                          if (value < -tolerance) {
    3087                               infeasible_->quickAdd(iSequence, value * value);
    3088                          }
    3089                     }
     3404                    infeas[iSequence] = 0.0;
     3405                    unsigned char thisStatus=status[iSequence]&7;
     3406                    if (thisStatus==3) {
     3407                    } else if ((thisStatus&1)!=0) {
     3408                      // basic or fixed
     3409                      value=0.0;
     3410                    } else if (thisStatus==2) {
     3411                      value=-value;
     3412                    } else {
     3413                      // free or superbasic
     3414                      if (fabs(value) > FREE_ACCEPT * tolerance) {
     3415                        // we are going to bias towards free (but only if reasonable)
     3416                        value = -fabs(value)*FREE_BIAS;
     3417                      } else {
     3418                        value=0.0;
     3419                      }
     3420                    }
     3421                    if (value < -tolerance) {
     3422                      infeas[iSequence] = value*value;
     3423                      index[numberNonZero++] = iSequence;
     3424                    }
    30903425               }
    30913426#else
     
    30933428               int numberColumns = model_->numberColumns();
    30943429               for (iSequence = 0; iSequence < numberColumns; iSequence++) {
     3430                    infeas[iSequence] = 0.0;
    30953431                    double value = reducedCost[iSequence];
    3096                     ClpSimplex::Status status = model_->getStatus(iSequence);
    3097 
    3098                     switch(status) {
    3099 
    3100                     case ClpSimplex::basic:
    3101                     case ClpSimplex::isFixed:
    3102                          break;
    3103                     case ClpSimplex::isFree:
    3104                     case ClpSimplex::superBasic:
    3105                          if (fabs(value) > FREE_ACCEPT * tolerance) {
    3106                               // check hasn't slipped through
    3107                               if (solution[iSequence]<lower[iSequence]+primalTolerance) {
    3108                                 model_->setStatus(iSequence,ClpSimplex::atLowerBound);
    3109                                 if (value < -tolerance) {
    3110                                   infeasible_->quickAdd(iSequence, value * value);
    3111                                 }
    3112                               } else if (solution[iSequence]>upper[iSequence]-primalTolerance) {
    3113                                 model_->setStatus(iSequence,ClpSimplex::atUpperBound);
    3114                                 if (value > tolerance) {
    3115                                   infeasible_->quickAdd(iSequence, value * value);
    3116                                 }
    3117                               } else {
    3118                                     // we are going to bias towards free (but only if reasonable)
    3119                                     value *= FREE_BIAS;
    3120                                     // store square in list
    3121                                     infeasible_->quickAdd(iSequence, value * value);
    3122                                   }
    3123                          }
    3124                          break;
    3125                     case ClpSimplex::atUpperBound:
    3126                          if (value > tolerance) {
    3127                               infeasible_->quickAdd(iSequence, value * value);
    3128                          }
    3129                          break;
    3130                     case ClpSimplex::atLowerBound:
    3131                          if (value < -tolerance) {
    3132                               infeasible_->quickAdd(iSequence, value * value);
    3133                          }
    3134                     }
     3432                    unsigned char thisStatus=status[iSequence]&7;
     3433                    if (thisStatus==3) {
     3434                    } else if ((thisStatus&1)!=0) {
     3435                      // basic or fixed
     3436                      value=0.0;
     3437                    } else if (thisStatus==2) {
     3438                      value=-value;
     3439                    } else {
     3440                      // free or superbasic
     3441                      if (fabs(value) > FREE_ACCEPT * tolerance) {
     3442                        // check hasn't slipped through
     3443                        if (solution[iSequence]<lower[iSequence]+primalTolerance) {
     3444                          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     3445                        } else if (solution[iSequence]>upper[iSequence]-primalTolerance) {
     3446                          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
     3447                          value = - value;
     3448                        } else {
     3449                          // we are going to bias towards free (but only if reasonable)
     3450                          value = -fabs(value)*FREE_BIAS;
     3451                        }
     3452                      } else {
     3453                        value=0.0;
     3454                      }
     3455                    }
     3456                    if (value < -tolerance) {
     3457                      infeas[iSequence] = value*value;
     3458                      index[numberNonZero++] = iSequence;
     3459                    }
    31353460               }
    31363461               // Rows
    31373462               for ( ; iSequence < number; iSequence++) {
    31383463                    double value = reducedCost[iSequence];
    3139                     ClpSimplex::Status status = model_->getStatus(iSequence);
    3140 
    3141                     switch(status) {
    3142 
    3143                     case ClpSimplex::basic:
    3144                     case ClpSimplex::isFixed:
    3145                          break;
    3146                     case ClpSimplex::isFree:
    3147                     case ClpSimplex::superBasic:
    3148                          if (fabs(value) > FREE_ACCEPT * tolerance) {
    3149                               // we are going to bias towards free (but only if reasonable)
    3150                               value *= FREE_BIAS;
    3151                               // store square in list
    3152                               infeasible_->quickAdd(iSequence, value * value);
    3153                          }
    3154                          break;
    3155                     case ClpSimplex::atUpperBound:
    3156                          if (value > tolerance) {
    3157                               infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
    3158                          }
    3159                          break;
    3160                     case ClpSimplex::atLowerBound:
    3161                          if (value < -tolerance) {
    3162                               infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
    3163                          }
    3164                     }
    3165                }
    3166 #endif
     3464                    infeas[iSequence]=0.0;
     3465                    unsigned char thisStatus=status[iSequence]&7;
     3466                    if (thisStatus==3) {
     3467                    } else if ((thisStatus&1)!=0) {
     3468                      // basic or fixed
     3469                      value=0.0;
     3470                    } else if (thisStatus==2) {
     3471                      value=-value;
     3472                    } else {
     3473                      // free or superbasic
     3474                      if (fabs(value) > FREE_ACCEPT * tolerance) {
     3475                        // we are going to bias towards free (but only if reasonable)
     3476                        value = -fabs(value)*FREE_BIAS;
     3477                      } else {
     3478                        value=0.0;
     3479                      }
     3480                    }
     3481                    if (value < -tolerance) {
     3482                      infeas[iSequence] = value*value;
     3483                      index[numberNonZero++] = iSequence;
     3484                    }
     3485               }
     3486#endif
     3487               infeasible_->setNumElements(numberNonZero);
    31673488          } else {
    31683489               ClpNonLinearCost * nonLinear = model_->nonLinearCost();
     3490               infeasible_->clear();
    31693491               // can go both ways
    31703492               for (iSequence = 0; iSequence < number; iSequence++) {
     
    32553577     int * which = input->getIndices();
    32563578     double * work = input->denseVector();
     3579#if ALT_UPDATE_WEIGHTS !=2
    32573580     int newNumber = 0;
    32583581     int * newWhich = alternateWeights_->getIndices();
    32593582     double * newWork = alternateWeights_->denseVector();
     3583#endif
     3584#ifdef ALT_UPDATE_WEIGHTSz
     3585     {
     3586       //int newNumber2 = 0;
     3587       if (!altVector[0]) {
     3588         altVector[0]=new CoinIndexedVector(2000);
     3589         altVector[1]=new CoinIndexedVector(2000);
     3590         altVector[2]=new CoinIndexedVector(2000);
     3591       }
     3592       altVector[0]->clear();
     3593       // get updated pivot row
     3594       int pivotRow = model_->pivotRow();
     3595       // should it be - or what
     3596       altVector[0]->quickAdd(pivotRow,model_->dualIn());
     3597       model_->factorization()->updateColumnTranspose(altVector[2],
     3598                                                      altVector[0]);
     3599       double * work2 = altVector[0]->denseVector();
     3600       //altVector[1]->clear();
     3601       int * newWhich2 = altVector[1]->getIndices();
     3602       double * newWork2 = altVector[1]->denseVector();
     3603       int number2 = altVector[1]->getNumElements();
     3604       int nRow=model_->numberRows();
     3605       int nCol=model_->numberColumns();
     3606       int nTotal=nRow+nCol;
     3607       double * temp=new double[2*nTotal+nRow];
     3608       memset(temp,0,(2*nTotal+nRow)*sizeof(double));
     3609       double * pivRow = temp+nTotal;
     3610       double * temp2 = temp+nCol;
     3611       double * temp2P=pivRow+nCol;
     3612       double * piU=pivRow+nTotal;
     3613       double devex=0.0;
     3614       double scaleFactor = 1.0/model_->dualIn();
     3615       const int * pivotVariable = model_->pivotVariable();
     3616       for (int i = 0; i < number; i++) {
     3617         int iRow = which[i];
     3618         int iPivot = pivotVariable[iRow];
     3619         if (reference(iPivot)) {
     3620           devex += work[iRow] * work[iRow];
     3621         }
     3622       }
     3623       int sequenceIn = model_->sequenceIn();
     3624       int sequenceOut = model_->sequenceOut();
     3625       for (int i=0;i<number2;i++) {
     3626         int iRow=newWhich2[i];
     3627         temp2[iRow] = newWork2[iRow];
     3628       }
     3629       //if (!input->packedMode()) {
     3630       for (int i=0;i<number;i++) {
     3631         int iRow=which[i];
     3632         piU[iRow] = work2[iRow];
     3633         temp2P[iRow]=work2[iRow];
     3634       }
     3635       double alpha=model_->alpha();
     3636       const int * row = model_->matrix()->getIndices();
     3637       const CoinBigIndex * columnStart = model_->matrix()->getVectorStarts();
     3638       const int * columnLength = model_->matrix()->getVectorLengths();
     3639       const double * element = model_->matrix()->getElements();
     3640       for (int i=0;i<nCol;i++) {
     3641         CoinBigIndex start = columnStart[i];
     3642         CoinBigIndex end = start + columnLength[i];
     3643         double value=0.0;
     3644         double value2=0.0;
     3645         for (CoinBigIndex j = start; j < end; j++) {
     3646           int iRow = row[j];
     3647           value -= piU[iRow] * element[j];
     3648           value2 -= newWork2[iRow] * element[j];
     3649         }
     3650         pivRow[i]=value;
     3651         temp[i]=value2;
     3652       }
     3653       const unsigned char * COIN_RESTRICT statusArray = model_->statusArray();
     3654       for (int i=0;i<nTotal;i++) {
     3655        unsigned char thisStatus=statusArray[i]&7;
     3656#if 0
     3657        if (thisStatus==3) {
     3658          // lower
     3659        } else if ((thisStatus&1)!=0) {
     3660          // basic or fixed
     3661          //value=0.0;
     3662        } else if (thisStatus==2) {
     3663          // uppervalue=-value;
     3664        } else {
     3665          // free or superbasic
     3666          //if (fabs(value) > FREE_ACCEPT * -dualTolerance) {
     3667            // we are going to bias towards free (but only if reasonable)
     3668            //value = -fabs(value)*FREE_BIAS;
     3669          //} else {
     3670          //value=0.0;
     3671          //}
     3672        }
     3673#else
     3674        if (thisStatus!=1&&thisStatus!=5) {
     3675          double pivot = pivRow[i]*scaleFactor;
     3676          double modification = temp[i];
     3677          double thisWeight=weights_[i];
     3678          double pivotSquared = pivot*pivot;
     3679          double newWeight = thisWeight + pivotSquared * devexA - 2.0*pivot * modification;
     3680          temp[i]=newWeight;
     3681        } else {
     3682          temp[i]=1.0;
     3683        }
     3684#endif
     3685       }
     3686       temp[sequenceOut]=devexA/(alpha*alpha);
     3687       // to keep clean for debug
     3688#ifndef NDEBUG
     3689       {
     3690         if (sequenceOut<nCol) {
     3691           if (model_->columnLower()[sequenceOut]==
     3692               model_->columnUpper()[sequenceOut])
     3693             temp[sequenceOut]=1.0;
     3694         } else {
     3695           if (model_->rowLower()[sequenceOut-nCol]==
     3696               model_->rowUpper()[sequenceOut-nCol])
     3697             temp[sequenceOut]=1.0;
     3698         }
     3699       }
     3700#endif
     3701       temp[sequenceIn]=1.0;
     3702       for (int i=0;i<nTotal;i++) {
     3703         printf("%g ",temp[i]);
     3704         if ((i%10)==9)
     3705           printf("\n");
     3706       }
     3707       if (((nTotal-1)%10)!=9)
     3708         printf("\n");
     3709       delete [] temp;
     3710     }
     3711#endif
    32603712     int i;
    32613713     int sequenceIn = model_->sequenceIn();
     
    32743726                         int iRow = which[i];
    32753727                         devex_ += work[iRow] * work[iRow];
     3728#if ALT_UPDATE_WEIGHTS !=2
    32763729                         newWork[iRow] = -2.0 * work[iRow];
    3277                     }
     3730#endif
     3731                    }
     3732#if ALT_UPDATE_WEIGHTS !=2
    32783733                    newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
     3734#endif
    32793735                    devex_ += ADD_ONE;
    32803736                    weights_[sequenceOut] = 1.0 + ADD_ONE;
     3737#if ALT_UPDATE_WEIGHTS !=2
    32813738                    CoinMemcpyN(which, number, newWhich);
    32823739                    alternateWeights_->setNumElements(number);
     3740#endif
    32833741               } else {
    32843742                    if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
     
    32883746                              if (reference(iPivot)) {
    32893747                                   devex_ += work[iRow] * work[iRow];
     3748#if ALT_UPDATE_WEIGHTS !=2
    32903749                                   newWork[iRow] = -2.0 * work[iRow];
    32913750                                   newWhich[newNumber++] = iRow;
     3751#endif
    32923752                              }
    32933753                         }
     3754#if ALT_UPDATE_WEIGHTS !=2
    32943755                         if (!newWork[pivotRow] && devex_ > 0.0)
    32953756                              newWhich[newNumber++] = pivotRow; // add if not already in
    32963757                         newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
     3758#endif
    32973759                    } else {
    32983760                         for (i = 0; i < number; i++) {
     
    33123774                         weights_[sequenceOut] = 1.0;
    33133775                    }
     3776#if ALT_UPDATE_WEIGHTS !=2
    33143777                    alternateWeights_->setNumElements(newNumber);
     3778#endif
    33153779               }
    33163780          } else {
     
    33383802               if (switchType == 1) {
    33393803                    for (i = 0; i < number; i++) {
     3804#if ALT_UPDATE_WEIGHTS !=2
    33403805                         int iRow = which[i];
     3806#endif
    33413807                         devex_ += work[i] * work[i];
     3808#if ALT_UPDATE_WEIGHTS !=2
    33423809                         newWork[iRow] = -2.0 * work[i];
    3343                     }
     3810#endif
     3811                    }
     3812#if ALT_UPDATE_WEIGHTS !=2
    33443813                    newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
     3814#endif
    33453815                    devex_ += ADD_ONE;
    33463816                    weights_[sequenceOut] = 1.0 + ADD_ONE;
     3817#if ALT_UPDATE_WEIGHTS !=2
    33473818                    CoinMemcpyN(which, number, newWhich);
    33483819                    alternateWeights_->setNumElements(number);
     3820#endif
    33493821               } else {
    33503822                    if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
     
    33543826                              if (reference(iPivot)) {
    33553827                                   devex_ += work[i] * work[i];
     3828#if ALT_UPDATE_WEIGHTS !=2
    33563829                                   newWork[iRow] = -2.0 * work[i];
    33573830                                   newWhich[newNumber++] = iRow;
     3831#endif
    33583832                              }
    33593833                         }
     3834#if ALT_UPDATE_WEIGHTS !=2
    33603835                         if (!newWork[pivotRow] && devex_ > 0.0)
    33613836                              newWhich[newNumber++] = pivotRow; // add if not already in
    33623837                         newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
     3838#endif
    33633839                    } else {
    33643840                         for (i = 0; i < number; i++) {
     
    33783854                         weights_[sequenceOut] = 1.0;
    33793855                    }
     3856#if ALT_UPDATE_WEIGHTS !=2
    33803857                    alternateWeights_->setNumElements(newNumber);
     3858#endif
    33813859               }
    33823860          } else {
     
    34283906                         << CoinMessageEol;
    34293907               initializeWeights();
     3908               // redo devex_
     3909               if (pivotRow >= 0)
     3910                 devex_=1.0;
    34303911          }
    34313912     }
     
    34743955     }
    34753956
    3476      double oldDevex = weights_[sequence];
     3957     double oldDevex = CoinMax(weights_[sequence],1.0e-4);
     3958     devex = CoinMax(devex,1.0e-4);
    34773959     double check = CoinMax(devex, oldDevex);;
     3960     rowArray1->setNumElements(0);
    34783961     if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
    3479        COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
     3962       //COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
     3963       printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex);
     3964       if (mode_==0) {
     3965         rowArray1->setNumElements(0);
     3966         model_->unpack(rowArray1, sequence);
     3967         number = rowArray1->getNumElements();
     3968         for (i = 0; i< number;i++)
     3969           printf("(%d,%g) ",which[i],work[which[i]]);
     3970         printf("\n");
     3971         model_->factorization()->updateColumn(rowArray2, rowArray1);
     3972         number = rowArray1->getNumElements();
     3973         for (i = 0; i< number;i++)
     3974           printf("(%d,%g) ",which[i],work[which[i]]);
     3975         printf("\n");
     3976         devex=0.0;
     3977         for (i = 0; i < number; i++) {
     3978           int iRow = which[i];
     3979           int iPivot = pivotVariable[iRow];
     3980           if (reference(iPivot)) {
     3981             devex += work[iRow] * work[iRow];
     3982           }
     3983           work[iRow] = 0.0;
     3984         }
     3985         if (reference(sequence))
     3986           devex += 1.0;
     3987       }
    34803988          // update so won't print again
    34813989          weights_[sequence] = devex;
    34823990     }
    3483      rowArray1->setNumElements(0);
    34843991}
    34853992
  • trunk/Clp/src/ClpPrimalColumnSteepest.hpp

    r2149 r2235  
    6161                      CoinIndexedVector * spareColumn1,
    6262                      CoinIndexedVector * spareColumn2);
    63      /// Update djs, weights for Steepest using djs
     63     /** Update djs, weights for Steepest using djs
     64         sets best sequence (possibly) */
    6465     void djsAndSteepest(CoinIndexedVector * updates,
    6566                         CoinIndexedVector * spareRow2,
     
    8788                       CoinIndexedVector * spareColumn2);
    8889     /// Updates two arrays for steepest
    89      void transposeTimes2(const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
     90     int transposeTimes2(const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    9091                          const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
    9192                          CoinIndexedVector * spare, double scaleFactor);