Changeset 2078 for trunk


Ignore:
Timestamp:
Jan 5, 2015 7:39:49 AM (5 years ago)
Author:
forrest
Message:

make faster on difficult problems

Location:
trunk/Clp/src
Files:
27 edited

Legend:

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

    r2025 r2078  
    8282          }
    8383          if (rhs.alternateWeights_) {
     84#ifndef ABC_USE_COIN_FACTORIZATIONz
    8485               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     86               // zero
     87               CoinZeroN(alternateWeights_->getIndices(),
     88                         alternateWeights_->capacity());
     89#else
     90#endif
    8591          } else {
    8692               alternateWeights_ = NULL;
     
    152158          }
    153159          if (rhs.alternateWeights_ != NULL) {
     160#ifndef ABC_USE_COIN_FACTORIZATIONz
    154161               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     162               // zero
     163               CoinZeroN(alternateWeights_->getIndices(),
     164                         alternateWeights_->capacity());
     165#else
     166#endif
    155167          } else {
    156168               alternateWeights_ = NULL;
     
    13951407          alternateWeights_ = new CoinIndexedVector();
    13961408          // enough space so can use it for factorization
    1397           // enoughfor ordered
     1409#ifndef ABC_USE_COIN_FACTORIZATION
     1410          // enough for ordered
    13981411          alternateWeights_->reserve(2*model_->numberRows() +
    13991412                                     model_->factorization()->maximumPivots());
     1413#else
     1414          int n=(2*model_->numberRows()+
     1415             model_->factorization()->maximumPivots()+7)&~3;
     1416          alternateWeights_->reserve(n);
     1417          // zero
     1418          CoinZeroN(alternateWeights_->getIndices(),
     1419                    alternateWeights_->capacity());
     1420#endif
    14001421     }
    14011422}
     
    14661487        alternateWeights_ = new CoinIndexedVector();
    14671488        // enough space so can use it for factorization
     1489#ifndef ABC_USE_COIN_FACTORIZATION
    14681490        alternateWeights_->reserve(2*numberRows +
    14691491                                   model_->factorization()->maximumPivots());
     1492#else
     1493        int n=(2*model_->numberRows()+
     1494           model_->factorization()->maximumPivots()+7)&~3;
     1495        alternateWeights_->reserve(n);
     1496        // zero
     1497        CoinZeroN(alternateWeights_->getIndices(),
     1498                  alternateWeights_->capacity());
     1499#endif
    14701500        initializeWeights();
    14711501        // create saved weights
     
    19031933     } else {
    19041934          CoinIndexedVector * temp = new CoinIndexedVector();
     1935#ifndef ABC_USE_COIN_FACTORIZATIONz
    19051936          temp->reserve(numberRows +
    19061937                        model_->factorization()->maximumPivots());
     1938#else
     1939#endif
    19071940          double * array = alternateWeights_->denseVector();
    19081941          int * which = alternateWeights_->getIndices();
  • trunk/Clp/src/AbcSimplex.cpp

    r2074 r2078  
    12551255    }
    12561256  }
     1257#ifdef ABC_USE_COIN_FACTORIZATION
     1258  // sparse methods
     1259  abcFactorization_->goSparse();
     1260#endif
    12571261  return status;
    12581262}
     
    15001504    }
    15011505    //printf("ZZ maxpivots %d its %d\n",numberPivots,maximumPivots);
     1506#if CLP_FACTORIZATION_NEW_TIMING>1
     1507    abcFactorization_->statsRefactor('M');
     1508#endif
    15021509    return 1;
    15031510  } else if ((abcFactorization_->timeToRefactorize() && !dontInvert)
    15041511             ||invertNow) {
    15051512    //printf("ZZ time invertNow %s its %d\n",invertNow ? "yes":"no",numberPivots);
     1513#if CLP_FACTORIZATION_NEW_TIMING>1
     1514    abcFactorization_->statsRefactor('T');
     1515#endif
    15061516    return 1;
    15071517  } else if (forceFactorization_ > 0 &&
     
    15171527      forceFactorization_ = -1; //off
    15181528    //printf("ZZ forceFactor %d its %d\n",forceFactorization_,numberPivots);
     1529#if CLP_FACTORIZATION_NEW_TIMING>1
     1530    abcFactorization_->statsRefactor('F');
     1531#endif
    15191532    return 1;
    15201533  } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
  • trunk/Clp/src/AbcSimplexDual.cpp

    r2070 r2078  
    9797#include "ClpEventHandler.hpp"
    9898#include "AbcSimplexFactorization.hpp"
     99#include "AbcNonLinearCost.hpp"
    99100#include "CoinPackedMatrix.hpp"
    100101#include "CoinIndexedVector.hpp"
     
    23372338AbcSimplex::gutsOfSolution(int type)
    23382339{
     2340  double oldLargestPrimalError=largestPrimalError_;
     2341  double oldLargestDualError=largestDualError_;
    23392342  AbcSimplexDual * dual = reinterpret_cast<AbcSimplexDual *>(this);
    23402343  // do work and check
     
    24212424    //if (abcFactorization_->zeroTolerance() > 1.0e-18)
    24222425    //abcFactorization_->zeroTolerance(1.0e-18);
     2426  }
     2427  bool notChanged=true;
     2428  if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_<0 || abcFactorization_->pivotTolerance()<0.9899999999) &&
     2429      (oldLargestDualError||oldLargestPrimalError)) {
     2430    double useOldDualError=oldLargestDualError;
     2431    double useDualError=largestDualError_;
     2432    if (algorithm_>0&&abcNonLinearCost_&&
     2433        abcNonLinearCost_->sumInfeasibilities()) {
     2434      double factor=CoinMax(1.0,CoinMin(1.0e3,infeasibilityCost_*1.0e-6));
     2435      useOldDualError /= factor;
     2436      useDualError /= factor;
     2437    }
     2438    if ((largestPrimalError_>1.0e3&&
     2439         oldLargestPrimalError*1.0e2<largestPrimalError_)||
     2440        (useDualError>1.0e3&&
     2441         useOldDualError*1.0e2<useDualError)) {
     2442      double pivotTolerance = abcFactorization_->pivotTolerance();
     2443      double factor=(largestPrimalError_>1.0e10||largestDualError_>1.0e10)
     2444        ? 2.0 : 1.2;
     2445      if (pivotTolerance<0.1)
     2446        abcFactorization_->pivotTolerance(0.1);
     2447      else if (pivotTolerance<0.98999999)
     2448        abcFactorization_->pivotTolerance(CoinMin(0.99,pivotTolerance*factor));
     2449      notChanged=pivotTolerance==abcFactorization_->pivotTolerance();
     2450#if defined(CLP_USEFUL_PRINTOUT) && !defined(GCC_4_9)
     2451      if (pivotTolerance<0.9899999) {
     2452        double newTolerance=abcFactorization_->pivotTolerance();
     2453        printf("Changing pivot tolerance from %g to %g and backtracking\n",
     2454               pivotTolerance,newTolerance);
     2455      }
     2456      printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
     2457             oldLargestPrimalError,largestPrimalError_,
     2458             oldLargestDualError,largestDualError_,
     2459             pivotTolerance);
     2460#endif
     2461      if (pivotTolerance<0.9899999) {
     2462        //largestPrimalError_=0.0;
     2463        //largestDualError_=0.0;
     2464        //returnCode=1;
     2465      }
     2466    }
     2467  }
     2468  if (progress_.iterationNumber_[0]>0&&
     2469      progress_.iterationNumber_[CLP_PROGRESS-1]
     2470      -progress_.iterationNumber_[0]<CLP_PROGRESS*3&&
     2471      abcFactorization_->pivotTolerance()<0.25&&notChanged) {
     2472    double pivotTolerance = abcFactorization_->pivotTolerance();
     2473    abcFactorization_->pivotTolerance(pivotTolerance*1.5);
     2474#if defined(CLP_USEFUL_PRINTOUT) && !defined(GCC_4_9)
     2475    double newTolerance=abcFactorization_->pivotTolerance();
     2476    printf("Changing pivot tolerance from %g to %g - inverting too often\n",
     2477           pivotTolerance,newTolerance);
     2478#endif
    24232479  }
    24242480  return numberRefinements;
  • trunk/Clp/src/AbcSimplexFactorization.cpp

    r1910 r2078  
    1717#include "AbcMatrix.hpp"
    1818#include "CoinAbcFactorization.hpp"
    19 //#include "CoinDenseFactorization.hpp"
     19#include "CoinFactorization.hpp"
    2020#ifdef ABC_JUST_ONE_FACTORIZATION
    2121#define CoinAbcFactorization CoinAbcBaseFactorization
     
    3131#undef CoinAbcSmallFactorization
    3232#define CoinAbcSmallFactorization CoinAbcOrderedFactorization
     33#endif
     34#ifdef ABC_USE_COIN_FACTORIZATION
     35#undef CoinAbcFactorization
     36#undef CoinAbcSmallFactorization
     37#undef CoinAbcLongFactorization
     38#undef CoinAbcOrderedFactorization
     39#define CoinAbcFactorization CoinFactorization
     40#define CoinAbcSmallFactorization CoinFactorization
     41#define CoinAbcLongFactorization CoinFactorization
     42#define CoinAbcOrderedFactorization CoinFactorization
     43#endif
     44#ifndef ABC_USE_COIN_FACTORIZATION
     45#undef CLP_FACTORIZATION_NEW_TIMING
     46#else
     47#ifndef CLP_FACTORIZATION_NEW_TIMING
     48#define CLP_FACTORIZATION_NEW_TIMING 1
     49#endif
    3350#endif
    3451//-------------------------------------------------------------------
     
    5774  goLongThreshold_ = rhs.goLongThreshold_;
    5875  numberSlacks_=rhs.numberSlacks_;
     76  model_=rhs.model_;
     77#ifndef ABC_USE_COIN_FACTORIZATION
    5978  int goDense = 0;
    60   model_=rhs.model_;
    6179  if (denseIfSmaller > 0 && denseIfSmaller <= goDenseThreshold_) {
    6280    CoinAbcDenseFactorization * denseR =
     
    99117    coinAbcFactorization_->zeroTolerance(rhs.coinAbcFactorization_->zeroTolerance());
    100118  }
     119#else
     120  coinAbcFactorization_ = new CoinFactorization(*rhs.coinAbcFactorization_);
     121#endif
    101122}
    102123//-------------------------------------------------------------------
     
    124145    if (rhs.coinAbcFactorization_) {
    125146      delete coinAbcFactorization_;
     147#ifndef ABC_USE_COIN_FACTORIZATION
    126148      coinAbcFactorization_ = rhs.coinAbcFactorization_->clone();
     149#else
     150      coinAbcFactorization_ = new CoinFactorization(*rhs.coinAbcFactorization_);
     151#endif
    127152    }
    128153  }
     
    133158AbcSimplexFactorization::goDenseOrSmall(int numberRows)
    134159{
     160#ifndef ABC_USE_COIN_FACTORIZATION
    135161  if (!forceB_) {
    136162    delete coinAbcFactorization_;
     
    145171    }
    146172  }
     173#endif
    147174}
    148175// If nonzero force use of 1,dense 2,small 3,long
     
    150177AbcSimplexFactorization::forceOtherFactorization(int which)
    151178{
     179#ifndef ABC_USE_COIN_FACTORIZATION
    152180  delete coinAbcFactorization_;
    153181  forceB_ = 0;
     
    174202    coinAbcFactorization_ = new CoinAbcFactorization();
    175203  }
     204#endif
     205}
     206#ifdef CLP_FACTORIZATION_NEW_TIMING
     207static bool readTwiddle=false;
     208static double weightIncU=1.0;
     209static double weightR=2.0;
     210static double weightRest=1.0;
     211static double weightFactL=30.0;
     212static double weightFactDense=0.1;
     213static double weightNrows=10.0;
     214static double increaseNeeded=1.1;
     215static double constWeightIterate = 1.0;
     216static double weightNrowsIterate = 3.0;
     217bool
     218AbcSimplexFactorization::timeToRefactorize() const
     219{
     220  bool reFactor = (coinAbcFactorization_->pivots() * 3 > coinAbcFactorization_->maximumPivots() * 2 &&
     221                   coinAbcFactorization_->numberElementsR() * 3 > (coinAbcFactorization_->numberElementsL() +
     222                                                                   coinAbcFactorization_->numberElementsU()) * 2 + 1000 &&
     223                   !coinAbcFactorization_->numberDense());
     224  bool reFactor3=false;
     225  int numberPivots=coinAbcFactorization_->pivots();
     226  //if (coinAbcFactorization_->pivots()<2)
     227  if (numberPivots>lastNumberPivots_) {
     228    if (!lastNumberPivots_) {
     229      //lastR=0;
     230      //lastU=endLengthU;
     231      totalInR_=0.0;
     232      totalInIncreasingU_=0.0;
     233      shortestAverage_=COIN_DBL_MAX;
     234      if (!readTwiddle) {
     235        readTwiddle=true;
     236        char * environ = getenv("CLP_TWIDDLE");
     237        if (environ) {
     238          sscanf(environ,"%lg %lg %lg %lg %lg %lg %lg %lg %lg",
     239                 &weightIncU,&weightR,&weightRest,&weightFactL,
     240                 &weightFactDense,&weightNrows,&increaseNeeded,
     241                 &constWeightIterate,&weightNrowsIterate);
     242        }
     243        printf("weightIncU %g, weightR %g, weightRest %g, weightFactL %g, weightFactDense %g, weightNrows %g increaseNeeded %g constWeightIterate %g weightNrowsIterate %g\n",
     244               weightIncU,weightR,weightRest,weightFactL,
     245               weightFactDense,weightNrows,increaseNeeded,
     246               constWeightIterate,weightNrowsIterate);
     247      }
     248    }
     249    lastNumberPivots_=numberPivots;
     250    int numberDense=coinAbcFactorization_->numberDense();
     251    double nnd=numberDense*numberDense;
     252    int lengthL=coinAbcFactorization_->numberElementsL();
     253    int lengthR=coinAbcFactorization_->numberElementsR();
     254    int numberRows = coinAbcFactorization_->numberRows();
     255    int lengthU=coinAbcFactorization_->numberElementsU()-
     256      (numberRows-numberDense);
     257    totalInR_ += lengthR;
     258    int effectiveU=lengthU-effectiveStartNumberU_;
     259    totalInIncreasingU_ += effectiveU;
     260    //lastR=lengthR;
     261    //lastU=lengthU;
     262    double rest=lengthL+0.05*nnd;
     263    double constWeightFactor = weightFactL*lengthL+weightFactDense*nnd
     264      + weightNrows*numberRows;
     265    double constWeightIterateX = constWeightIterate*(lengthL+endLengthU_)
     266      + weightNrowsIterate*numberRows;
     267    double variableWeight = weightIncU*totalInIncreasingU_+
     268      weightR*totalInR_+weightRest*rest;
     269    double average=constWeightIterateX+
     270      (constWeightFactor+variableWeight)/static_cast<double>(numberPivots);
     271#if 0
     272    if ((numberPivots%20)==0&&!ifPrint3)
     273      printf("PIV %d nrow %d startU %d now %d L %d R %d dense %g average %g\n",
     274             numberPivots,numberRows,effectiveStartNumberU_,
     275             lengthU,lengthL,lengthR,nnd,average);
     276#endif
     277    shortestAverage_=CoinMin(shortestAverage_,average);
     278    if (average>increaseNeeded*shortestAverage_&&
     279        coinAbcFactorization_->pivots()>30) {
     280      //printf("PIVX %d nrow %d startU %d now %d L %d R %d dense %g average %g\n",
     281      //     numberPivots,numberRows,effectiveStartNumberU_,
     282      //     lengthU,lengthL,lengthR,nnd,average);
     283      reFactor3=true;
     284    }
     285  }
     286  if (reFactor|| reFactor3) {
     287    reFactor=true;
     288  }
     289  return reFactor;
     290}
     291#if CLP_FACTORIZATION_NEW_TIMING>1
     292void
     293AbcSimplexFactorization::statsRefactor(char when) const
     294{
     295  int numberPivots=coinAbcFactorization_->pivots();
     296  int numberDense=coinAbcFactorization_->numberDense();
     297  double nnd=numberDense*numberDense;
     298  int lengthL=coinAbcFactorization_->numberElementsL();
     299  int lengthR=coinAbcFactorization_->numberElementsR();
     300  int numberRows = coinAbcFactorization_->numberRows();
     301  int lengthU=coinAbcFactorization_->numberElementsU()-
     302    (numberRows-numberDense);
     303  double rest=lengthL+0.05*nnd;
     304  double constWeightFactor = weightFactL*lengthL+weightFactDense*nnd
     305    + weightNrows*numberRows;
     306  double constWeightIterateX = constWeightIterate*(lengthL+endLengthU_)
     307    + weightNrowsIterate*numberRows;
     308  double variableWeight = weightIncU*totalInIncreasingU_+
     309    weightR*totalInR_+weightRest*rest;
     310  double average=constWeightIterateX+
     311    (constWeightFactor+variableWeight)/static_cast<double>(numberPivots);
     312  printf("APIV%c %d nrow %d startU %d now %d L %d R %d dense %g average %g - shortest %g\n",
     313         when,numberPivots,numberRows,effectiveStartNumberU_,
     314         lengthU,lengthL,lengthR,nnd,average,shortestAverage_);
     315}
     316#endif
     317#else
     318bool
     319AbcSimplexFactorization::timeToRefactorize() const
     320{
     321  return coinAbcFactorization_->pivots() > coinAbcFactorization_->numberRows() / 2.45 + 20;
     322}
     323#endif
     324/* returns empty fake vector carved out of existing
     325   later - maybe use associated arrays */
     326static CoinIndexedVector * fakeVector(CoinIndexedVector * vector,
     327                                      int fakeCapacity)
     328{
     329  int oldCapacity=vector->capacity();
     330  CoinIndexedVector * newVector = new CoinIndexedVector();
     331  newVector->setCapacity(fakeCapacity);
     332  newVector->setDenseVector(vector->denseVector()+oldCapacity);
     333  newVector->setIndexVector(vector->getIndices()+
     334                            oldCapacity+((oldCapacity+3)>>2));
     335  vector->checkClean();
     336  newVector->checkClear();
     337  return newVector;
     338}
     339static void deleteFakeVector(CoinIndexedVector * vector,
     340                             CoinIndexedVector * fakeVector)
     341{
     342  int * index = vector->getIndices();
     343  fakeVector->checkClear();
     344  fakeVector->setCapacity(0);
     345  fakeVector->setDenseVector(NULL);
     346  fakeVector->setIndexVector(NULL);
     347  delete fakeVector;
     348  vector->checkClean();
    176349}
    177350// Synchronize stuff
     
    185358  goDenseOrSmall(model->numberRows());
    186359  maximumPivots(static_cast<int>(otherFactorization->maximumPivots()*1.2));
    187 }
     360#ifdef ABC_USE_COIN_FACTORIZATION
     361  // redo region sizes
     362  int maximumRows=model->numberRows()+maximumPivots()+1;
     363  int currentCapacity = model->usefulArray(0)->capacity();
     364  int newCapacity = currentCapacity+maximumRows+3;
     365  for (int i=0;i<ABC_NUMBER_USEFUL;i++) {
     366    CoinPartitionedVector * vector = model->usefulArray(i);
     367    vector->reserve(newCapacity);
     368    // zero
     369    CoinZeroN(vector->getIndices(),newCapacity);
     370    // but pretend old
     371    //vector->setCapacity(currentCapacity);
     372#if 0 //ndef NDEBUG
     373    vector->checkClear();
     374    CoinIndexedVector * newVector = fakeVector(vector,maximumRows);
     375    deleteFakeVector(vector,newVector);
     376#endif
     377  }
     378#endif
     379}
     380#ifdef CLP_FACTORIZATION_INSTRUMENT
     381extern double externalTimeStart;
     382extern double timeInFactorize;
     383extern double timeInUpdate;
     384extern double timeInUpdateTranspose;
     385extern double timeInUpdateFT;
     386extern double timeInUpdateTwoFT;
     387extern double timeInReplace;
     388extern int numberUpdate;
     389extern int numberUpdateTranspose;
     390extern int numberUpdateFT;
     391extern int numberUpdateTwoFT;
     392extern int numberReplace;
     393extern int currentLengthR;
     394extern int currentLengthU;
     395extern int currentTakeoutU;
     396#endif
    188397int
    189398AbcSimplexFactorization::factorize ( AbcSimplex * model,
    190399                                     int solveType, bool valuesPass)
    191400{
     401#ifdef CLP_FACTORIZATION_INSTRUMENT
     402  externalTimeStart=CoinCpuTime();
     403#endif
    192404  model_= model;
    193405  AbcMatrix * matrix = model->abcMatrix();
     
    197409  bool anyChanged = false;
    198410  coinAbcFactorization_->setStatus(-99);
    199   const int *  COIN_RESTRICT pivotVariable = model->pivotVariable();
     411#ifndef ABC_USE_COIN_FACTORIZATION
     412  const
     413#endif
     414    int *  COIN_RESTRICT pivotVariable = model->pivotVariable();
    200415  //returns 0 -okay, -1 singular, -2 too many in basis */
    201416  // allow dense
     
    241456    // Not needed for dense
    242457    numberElements = 3 * numberBasic + 3 * numberElements + 20000;
     458#ifndef ABC_USE_COIN_FACTORIZATION
    243459    int numberIterations = model->numberIterations();
    244460    coinAbcFactorization_->setUsefulInformation(&numberIterations, 0);
     461#else
     462    coinAbcFactorization_->gutsOfDestructor();
     463    coinAbcFactorization_->gutsOfInitialize(2);
     464#endif
    245465    coinAbcFactorization_->getAreas ( numberRows,
    246466                                    numberSlacks_ + numberColumnBasic, numberElements,
     
    257477    int *  COIN_RESTRICT numberInRow;
    258478    int *  COIN_RESTRICT numberInColumn;
     479#define slackValue 1.0
     480#ifndef ABC_USE_COIN_FACTORIZATION
    259481    elementU = coinAbcFactorization_->elements();
    260482    indexRowU = coinAbcFactorization_->indices();
    261483    startColumnU = coinAbcFactorization_->starts();
    262 #define slackValue 1.0
    263484    numberInRow = coinAbcFactorization_->numberInRow();
    264485    numberInColumn = coinAbcFactorization_->numberInColumn();
     
    266487    CoinZeroN ( numberInRow, numberRows  );
    267488    CoinZeroN ( numberInColumn, numberRows );
     489#else
     490    elementU = coinAbcFactorization_->elementU();
     491    indexRowU = coinAbcFactorization_->indexRowU();
     492    startColumnU = coinAbcFactorization_->startColumnU();
     493    numberInRow = coinAbcFactorization_->numberInRow();
     494    numberInColumn = coinAbcFactorization_->numberInColumn();
     495    CoinZeroN ( numberInRow, coinAbcFactorization_->numberRows() + 1 );
     496    CoinZeroN ( numberInColumn, coinAbcFactorization_->maximumColumnsExtra() + 1 );
     497#endif
    268498    for (i = 0; i < numberSlacks_; i++) {
    269499      int iRow = pivotTemp[i];
     
    286516    numberElements = startColumnU[numberRows-1]
    287517      + numberInColumn[numberRows-1];
     518#ifndef ABC_USE_COIN_FACTORIZATION
    288519    coinAbcFactorization_->preProcess ( );
    289520    coinAbcFactorization_->factor (model);
     521#else
     522    // recompute number basic
     523    numberBasic = numberSlacks_ + numberColumnBasic;
     524    if (numberBasic)
     525      numberElements = startColumnU[numberBasic-1]
     526        + numberInColumn[numberBasic-1];
     527    else
     528      numberElements = 0;
     529    coinAbcFactorization_->setNumberElementsU(numberElements);
     530#ifdef CLP_FACTORIZATION_NEW_TIMING
     531    lastNumberPivots_=0;
     532    effectiveStartNumberU_=numberElements-numberRows;
     533    //printf("%d slacks,%d in U at beginning\n",
     534    //numberRowBasic,numberElements);
     535#endif
     536    //saveFactorization("dump.d");
     537    if (coinAbcFactorization_->biasLU() >= 3 || coinAbcFactorization_->numberRows() != coinAbcFactorization_->numberColumns())
     538      coinAbcFactorization_->preProcess ( 2 );
     539    else
     540      coinAbcFactorization_->preProcess ( 3 ); // no row copy
     541    coinAbcFactorization_->factor (  );
     542#endif
    290543#if 0
    291544    if (model_->numberIterations()==23) {
     
    301554      coinAbcFactorization_->setSolveMode(solveMode);
    302555      coinAbcFactorization_->setStatus(-99);
    303     }
    304     if (coinAbcFactorization_->status() == -99)
     556    } else if (coinAbcFactorization_->status() == -99) {
     557      // get more memory
     558      coinAbcFactorization_->areaFactor( coinAbcFactorization_->areaFactor() * 2.0);
     559    }
     560    if (coinAbcFactorization_->status() == -99)
    305561      continue;
    306562    // If we get here status is 0 or -1
    307563    if (coinAbcFactorization_->status() == 0 && numberBasic == numberRows) {
     564#ifndef ABC_USE_COIN_FACTORIZATION
    308565      coinAbcFactorization_->postProcess(pivotTemp, model->pivotVariable());
     566#else
     567      const int * permuteBack = coinAbcFactorization_->permuteBack();
     568      const int * back = coinAbcFactorization_->pivotColumnBack();
     569      // Redo pivot order
     570      for (i = 0; i < numberRows; i++) {
     571        int k = pivotTemp[i];
     572        // so rowIsBasic[k] would be permuteBack[back[i]]
     573        int j = permuteBack[back[i]];
     574        //assert (pivotVariable[j] == -1);
     575        pivotVariable[j] = k;
     576      }
     577      // Set up permutation vector
     578      // these arrays start off as copies of permute
     579      // (and we could use permute_ instead of pivotColumn (not back though))
     580      ClpDisjointCopyN ( coinAbcFactorization_->permute(), numberRows , coinAbcFactorization_->pivotColumn()  );
     581      ClpDisjointCopyN ( coinAbcFactorization_->permuteBack(), numberRows , coinAbcFactorization_->pivotColumnBack()  );
     582      // See if worth going sparse and when
     583      coinAbcFactorization_->checkSparse();
     584#ifdef CLP_FACTORIZATION_NEW_TIMING
     585      endLengthU_ = coinAbcFactorization_->numberElements() -
     586        coinAbcFactorization_->numberDense()*coinAbcFactorization_->numberDense()
     587        -coinAbcFactorization_->numberElementsL();
     588#endif
     589#endif
    309590      model_->moveToBasic();
    310591    } else if (solveType == 0 || solveType == 2/*||solveType==1*/) {
     
    361642              } else {
    362643                thisStatus= AbcSimplex::isFixed;
    363                 solution[iPivot] = upper;
    364644              }
    365645            } else {
     
    425705  tempInfo[0] = model->numberIterations();
    426706  coinAbcFactorization_->setUsefulInformation(tempInfo, 1);
     707#ifdef CLP_FACTORIZATION_NEW_TIMING
     708  int nOld=0;
     709  int nNew=0;
     710  int seq;
     711  const CoinPackedMatrix * matrix=model->matrix();
     712  const int * columnLength = matrix->getVectorLengths();
     713  seq=model->sequenceIn();
     714  if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
     715    if (seq<model->numberRows())
     716      nNew=1;
     717    else
     718      nNew=columnLength[seq-model->numberRows()];
     719  }
     720  seq=model->sequenceOut();
     721  if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
     722    if (seq<model->numberRows())
     723      nOld=1;
     724    else
     725      nOld=columnLength[seq-model->numberRows()];
     726  }
     727  effectiveStartNumberU_ += nNew-nOld;
     728#endif
    427729  int returnCode =
    428730    coinAbcFactorization_->replaceColumn(tab ? tableauColumn : regionSparse,
     
    446748                                              double alpha )
    447749{
     750#ifndef ABC_USE_COIN_FACTORIZATION
    448751  bool tab = coinAbcFactorization_->wantsTableauColumn();
    449752  int tempInfo[1];
     
    458761                                              pivotRow,
    459762                                              alpha);
     763#else
     764#ifdef CLP_FACTORIZATION_NEW_TIMING
     765  int nOld=0;
     766  int nNew=0;
     767  int seq;
     768  const CoinPackedMatrix * matrix=model->matrix();
     769  const int * columnLength = matrix->getVectorLengths();
     770  seq=model->sequenceIn();
     771  if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
     772    if (seq<model->numberRows())
     773      nNew=1;
     774    else
     775      nNew=columnLength[seq-model->numberRows()];
     776  }
     777  seq=model->sequenceOut();
     778  if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
     779    if (seq<model->numberRows())
     780      nOld=1;
     781    else
     782      nOld=columnLength[seq-model->numberRows()];
     783  }
     784  effectiveStartNumberU_ += nNew-nOld;
     785#endif
     786  coinAbcFactorization_->replaceColumnPart3(regionSparse,
     787                                              pivotRow,
     788                                              alpha);
     789#endif
    460790}
    461791/* Replaces one Column to basis,
     
    472802                                              double alpha )
    473803{
     804#ifndef ABC_USE_COIN_FACTORIZATION
    474805  bool tab = coinAbcFactorization_->wantsTableauColumn();
    475806  int tempInfo[1];
     
    484815                                              pivotRow,
    485816                                              alpha);
     817#else
     818    coinAbcFactorization_->replaceColumnPart3(regionSparse,partialUpdate,
     819                                              pivotRow,
     820                                              alpha);
     821#endif
    486822}
    487823#if 0
     
    572908AbcSimplexFactorization::goSparse()
    573909{
     910#ifdef ABC_USE_COIN_FACTORIZATION
     911  // sparse methods
     912  coinAbcFactorization_->sparseThreshold(0);
     913  coinAbcFactorization_->goSparse();
     914#else
    574915  abort();
    575916  coinAbcFactorization_->goSparse();
     917#endif
    576918}
    577919// Set tolerances to safer of existing and given
  • trunk/Clp/src/AbcSimplexFactorization.hpp

    r2024 r2078  
    1515//#include "CoinAbcAnyFactorization.hpp"
    1616#include "AbcSimplex.hpp"
     17#ifndef ABC_USE_COIN_FACTORIZATION
    1718class ClpFactorization;
    18 
     19class CoinFactorization;
     20#else
     21#include "ClpFactorization.hpp"
     22#include "CoinFactorization.hpp"
     23#endif
    1924/** This just implements AbcFactorization when an AbcMatrix object
    2025    is passed.
     
    8893                                  int pivotRow)
    8994  {return coinAbcFactorization_->checkReplacePart1(regionSparse,partialUpdate,pivotRow);}
     95#ifdef MOVE_REPLACE_PART1A
    9096  /** Checks if can replace one Column to basis,
    9197      returns update alpha
     
    98104                             int pivotRow)
    99105  {return coinAbcFactorization_->checkReplacePart1b(regionSparse,pivotRow);}
     106#endif
    100107  /** Checks if can replace one Column to basis,
    101108      returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots */
     
    201208  /** Updates one column (FTRAN) */
    202209  inline void updateColumnCpu ( CoinIndexedVector & regionSparse,int whichCpu) const
     210#ifndef ABC_USE_COIN_FACTORIZATION
    203211  { coinAbcFactorization_->updateColumnCpu(regionSparse,whichCpu);}
     212#else
     213  { coinAbcFactorization_->updateColumn(regionSparse);}
     214#endif
    204215  /** Updates one column (BTRAN) */
    205216  inline void updateColumnTransposeCpu ( CoinIndexedVector & regionSparse,int whichCpu) const
     217#ifndef ABC_USE_COIN_FACTORIZATION
    206218  { coinAbcFactorization_->updateColumnTransposeCpu(regionSparse,whichCpu);}
     219#else
     220  { coinAbcFactorization_->updateColumnTranspose(regionSparse);}
     221#endif
    207222  /** Updates one full column (FTRAN) */
    208223  inline void updateFullColumn ( CoinIndexedVector & regionSparse) const
     
    213228  /** Updates one column for dual steepest edge weights (FTRAN) */
    214229  void updateWeights ( CoinIndexedVector & regionSparse) const
     230#ifndef ABC_USE_COIN_FACTORIZATION
    215231  { coinAbcFactorization_->updateWeights(regionSparse);}
     232#else
     233  { coinAbcFactorization_->updateColumn(regionSparse);}
     234#endif
    216235  //@}
    217236  /**@name Lifted from CoinFactorization */
     
    278297    return coinAbcFactorization_->numberDense() ;
    279298  }
    280   inline bool timeToRefactorize() const {
    281     return coinAbcFactorization_->pivots() > coinAbcFactorization_->numberRows() / 2.45 + 20;
    282   }
     299  bool timeToRefactorize() const;
     300#if CLP_FACTORIZATION_NEW_TIMING>1
     301  void statsRefactor(char when) const;
     302#endif
    283303  /// Get rid of all memory
    284304  inline void clearArrays() {
     
    362382  void goSparse();
    363383#ifndef NDEBUG
     384#ifndef ABC_USE_COIN_FACTORIZATION
    364385  inline void checkMarkArrays() const
    365386  { coinAbcFactorization_->checkMarkArrays();}
     387#else
     388  inline void checkMarkArrays() const
     389  { }
     390#endif
    366391#endif
    367392  /// Says whether to redo pivot order
    368393  inline bool needToReorder() const {abort();return true;}
    369394  /// Pointer to factorization
     395#ifndef ABC_USE_COIN_FACTORIZATION
    370396  CoinAbcAnyFactorization * factorization() const
    371397  { return coinAbcFactorization_;}
     398#else
     399  CoinFactorization * factorization() const
     400  { return coinAbcFactorization_;}
     401#endif
    372402  //@}
    373403 
     
    380410  AbcSimplex * model_;
    381411  /// Pointer to factorization
     412#ifndef ABC_USE_COIN_FACTORIZATION
    382413  CoinAbcAnyFactorization * coinAbcFactorization_;
     414#else
     415  CoinFactorization * coinAbcFactorization_;
     416#endif
     417#ifdef CLP_FACTORIZATION_NEW_TIMING
     418  /// For guessing when to re-factorize
     419  mutable double shortestAverage_;
     420  mutable double totalInR_;
     421  mutable double totalInIncreasingU_;
     422  mutable int endLengthU_;
     423  mutable int lastNumberPivots_;
     424  mutable int effectiveStartNumberU_;
     425#endif
    383426  /// If nonzero force use of 1,dense 2,small 3,long
    384427  int forceB_;
  • trunk/Clp/src/AbcSimplexPrimal.cpp

    r2074 r2078  
    298298    abcProgress_.fillFromModel(this);
    299299    abcProgress_.startCheck();
     300    bool needNewWeights=false;
    300301    double pivotTolerance=abcFactorization_->pivotTolerance();
    301302    /*
     
    331332      // may factorize, checks if problem finished
    332333      if (factorType)
    333         statusOfProblemInPrimal(factorType);
     334        statusOfProblemInPrimal(factorType + (needNewWeights ? 10 : 0));
     335      needNewWeights=false;
    334336      if (problemStatus_==10 && (moreSpecialOptions_ & 2048) != 0) {
    335337        problemStatus_=0;
     
    381383        break;
    382384#endif
    383       }
     385      } 
    384386     
    385387      // test for maximum iterations
     
    393395          // end of values pass
    394396          ifValuesPass = 0;
     397          needNewWeights=true;
    395398          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
    396399          if (status >= 0) {
     
    409412#endif
    410413        }
     414      }
     415      if (abcNonLinearCost_->numberInfeasibilities() > 0 && !abcProgress_.initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10
     416#ifdef TRY_SPLIT_VALUES_PASS
     417          && valuesStop<0
     418#endif
     419          ) {
     420        // first time infeasible - start up weight computation
     421        // compute with original costs
     422        CoinAbcMemcpy(djSaved_,abcDj_,numberTotal_);
     423        CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
     424        CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);
     425        gutsOfPrimalSolution(1);
     426        int numberSame = 0;
     427        int numberDifferent = 0;
     428        int numberZero = 0;
     429        int numberFreeSame = 0;
     430        int numberFreeDifferent = 0;
     431        int numberFreeZero = 0;
     432        int n = 0;
     433        for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
     434          if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
     435            // not basic
     436            double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
     437            double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
     438            double feasibleDj = abcDj_[iSequence];
     439            double infeasibleDj = djSaved_[iSequence] - feasibleDj;
     440            double value = feasibleDj * infeasibleDj;
     441            if (distanceUp > primalTolerance_) {
     442              // Check if "free"
     443              if (distanceDown > primalTolerance_) {
     444                // free
     445                if (value > dualTolerance_) {
     446                  numberFreeSame++;
     447                } else if(value < -dualTolerance_) {
     448                  numberFreeDifferent++;
     449                  abcDj_[n++] = feasibleDj / infeasibleDj;
     450                } else {
     451                  numberFreeZero++;
     452                }
     453              } else {
     454                // should not be negative
     455                if (value > dualTolerance_) {
     456                  numberSame++;
     457                } else if(value < -dualTolerance_) {
     458                  numberDifferent++;
     459                  abcDj_[n++] = feasibleDj / infeasibleDj;
     460                } else {
     461                  numberZero++;
     462                }
     463              }
     464            } else if (distanceDown > primalTolerance_) {
     465              // should not be positive
     466              if (value > dualTolerance_) {
     467                numberSame++;
     468              } else if(value < -dualTolerance_) {
     469                numberDifferent++;
     470                abcDj_[n++] = feasibleDj / infeasibleDj;
     471              } else {
     472                numberZero++;
     473              }
     474            }
     475      }
     476          abcProgress_.initialWeight_ = -1.0;
     477    }
     478#ifdef CLP_USEFUL_PRINTOUT
     479        printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
     480               numberSame,numberDifferent,numberZero,
     481               numberFreeSame,numberFreeDifferent,numberFreeZero);
     482#endif
     483        // we want most to be same
     484        if (n) {
     485          std::sort(abcDj_, abcDj_ + n);
     486#ifdef CLP_USEFUL_PRINTOUT
     487          double most = 0.95;
     488          int which = static_cast<int> ((1.0 - most) * static_cast<double> (n));
     489          double take2 = -abcDj_[which] * infeasibilityCost_;
     490          printf("XXXX inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take2,-abcDj_[0]*infeasibilityCost_,-abcDj_[n-1]*infeasibilityCost_);
     491#endif
     492          double take = -abcDj_[0] * infeasibilityCost_;
     493          // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);
     494          infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10);
     495#ifdef CLP_USEFUL_PRINTOUT
     496          printf("XXXX changing weight to %g\n",infeasibilityCost_);
     497#endif
     498          // may need to force redo of weights
     499          needNewWeights=true;
     500        }
     501        abcNonLinearCost_->checkInfeasibilities(0.0);
     502        gutsOfPrimalSolution(3);
    411503      }
    412504      // Check event
     
    721813#endif
    722814  //int saveType=type;
    723   if (type>=4) {
     815  if (type>3) {
    724816    // force factorization
    725     type &= 3;
    726817    numberPivots=9999999;
     818    if (type>10)
     819      type -=10;
     820    else
     821      type &= 3;
    727822  }
    728823#ifndef TRY_ABC_GUS
     
    821916    if (problemStatus_ != -4)
    822917      problemStatus_ = -3;
     918    if(abcProgress_.realInfeasibility_[0]<1.0e-1 &&
     919       primalTolerance_==1.0e-7&&abcProgress_.iterationNumber_[0]>0&&
     920       abcProgress_.iterationNumber_[CLP_PROGRESS-1]-abcProgress_.iterationNumber_[0]>25) {
     921      // default - so user did not set
     922      int iP;
     923      double minAverage=COIN_DBL_MAX;
     924      double maxAverage=0.0;
     925      for (iP=0;iP<CLP_PROGRESS;iP++) {
     926        int n=abcProgress_.numberInfeasibilities_[iP];
     927        if (!n) {
     928          break;
     929        } else {
     930          double average=abcProgress_.realInfeasibility_[iP];
     931          if (average>0.1)
     932            break;
     933          average /= static_cast<double>(n);
     934          minAverage=CoinMin(minAverage,average);
     935          maxAverage=CoinMax(maxAverage,average);
     936        }
     937      }
     938      if (iP==CLP_PROGRESS&&minAverage<1.0e-5&&maxAverage<1.0e-3) {
     939        // change tolerance
     940#if CBC_USEFUL_PRINTING>0
     941        printf("CCchanging tolerance\n");
     942#endif
     943        primalTolerance_=1.0e-6;
     944        dblParam_[ClpPrimalTolerance]=1.0e-6;
     945        moreSpecialOptions_ |= 4194304;
     946      }
     947    }
    823948    // at this stage status is -3 or -5 if looks unbounded
    824949    // get primal and dual solutions
     
    10911216    gutsOfPrimalSolution(2);
    10921217    abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
    1093   }
    1094   if (abcNonLinearCost_->numberInfeasibilities() > 0 && !abcProgress_.initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10
    1095 #ifdef TRY_SPLIT_VALUES_PASS
    1096       && valuesStop<0
    1097 #endif
    1098       ) {
    1099     // first time infeasible - start up weight computation
    1100     // compute with original costs
    1101     CoinAbcMemcpy(djSaved_,abcDj_,numberTotal_);
    1102     CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
    1103     gutsOfPrimalSolution(1);
    1104     int numberSame = 0;
    1105     int numberDifferent = 0;
    1106     int numberZero = 0;
    1107     int numberFreeSame = 0;
    1108     int numberFreeDifferent = 0;
    1109     int numberFreeZero = 0;
    1110     int n = 0;
    1111     for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
    1112       if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
    1113         // not basic
    1114         double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
    1115         double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
    1116         double feasibleDj = abcDj_[iSequence];
    1117         double infeasibleDj = djSaved_[iSequence] - feasibleDj;
    1118         double value = feasibleDj * infeasibleDj;
    1119         if (distanceUp > primalTolerance_) {
    1120           // Check if "free"
    1121           if (distanceDown > primalTolerance_) {
    1122             // free
    1123             if (value > dualTolerance_) {
    1124               numberFreeSame++;
    1125             } else if(value < -dualTolerance_) {
    1126               numberFreeDifferent++;
    1127               abcDj_[n++] = feasibleDj / infeasibleDj;
    1128             } else {
    1129               numberFreeZero++;
    1130             }
    1131           } else {
    1132             // should not be negative
    1133             if (value > dualTolerance_) {
    1134               numberSame++;
    1135             } else if(value < -dualTolerance_) {
    1136               numberDifferent++;
    1137               abcDj_[n++] = feasibleDj / infeasibleDj;
    1138             } else {
    1139               numberZero++;
    1140             }
    1141           }
    1142         } else if (distanceDown > primalTolerance_) {
    1143           // should not be positive
    1144           if (value > dualTolerance_) {
    1145             numberSame++;
    1146           } else if(value < -dualTolerance_) {
    1147             numberDifferent++;
    1148             abcDj_[n++] = feasibleDj / infeasibleDj;
    1149           } else {
    1150             numberZero++;
    1151           }
    1152         }
    1153       }
    1154       abcProgress_.initialWeight_ = -1.0;
    1155     }
    1156     //printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
    1157     //   numberSame,numberDifferent,numberZero,
    1158     //   numberFreeSame,numberFreeDifferent,numberFreeZero);
    1159     // we want most to be same
    1160     if (n) {
    1161       double most = 0.95;
    1162       std::sort(abcDj_, abcDj_ + n);
    1163       int which = static_cast<int> ((1.0 - most) * static_cast<double> (n));
    1164       double take = -abcDj_[which] * infeasibilityCost_;
    1165       //printf("XXXXZ inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take,-abcDj_[0]*infeasibilityCost_,-abcDj_[n-1]*infeasibilityCost_);
    1166       take = -abcDj_[0] * infeasibilityCost_;
    1167       infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);;
    1168       //printf("XXXX increasing weight to %g\n",infeasibilityCost_);
    1169     }
    1170     abcNonLinearCost_->checkInfeasibilities(0.0);
    1171     gutsOfPrimalSolution(3);
    11721218  }
    11731219  double trueInfeasibility = abcNonLinearCost_->sumInfeasibilities();
     
    41724218AbcSimplexPrimal::updatePartialUpdate(CoinIndexedVector & partialUpdate)
    41734219{
     4220#ifndef ABC_USE_COIN_FACTORIZATION
    41744221  CoinAbcFactorization * factorization=
    41754222    dynamic_cast<CoinAbcFactorization *>(abcFactorization_->factorization());
    41764223  if (factorization)
    41774224    factorization->updatePartialUpdate(partialUpdate);
     4225#else
     4226  abcFactorization_->factorization()->updatePartialUpdate(partialUpdate);
     4227#endif
    41784228}
    41794229// Create primal ray
  • trunk/Clp/src/CbcOrClpParam.cpp

    r2074 r2078  
    17901790     parameters[numberParameters-1].append("dw");
    17911791     parameters[numberParameters-1].append("idiot");
     1792#else
     1793     parameters[numberParameters-1].append("idiot1");
     1794     parameters[numberParameters-1].append("idiot2");
     1795     parameters[numberParameters-1].append("idiot3");
     1796     parameters[numberParameters-1].append("idiot4");
     1797     parameters[numberParameters-1].append("idiot5");
     1798     parameters[numberParameters-1].append("idiot6");
     1799     parameters[numberParameters-1].append("idiot7");
    17921800#endif
    17931801     parameters[numberParameters-1].setLonghelp
     
    30113019     parameters[numberParameters-1].setLonghelp
    30123020     (
    3013           "Normally Presolve does 5 passes but you may want to do less to make it\
     3021          "Normally Presolve does 10 passes but you may want to do less to make it\
    30143022 more lightweight or do more if improvements are still being made.  As Presolve will return\
    30153023 if nothing is being taken out, you should not normally need to use this fine tuning."
  • trunk/Clp/src/ClpFactorization.cpp

    r1903 r2078  
    2424#endif
    2525#endif
    26 //#define CLP_FACTORIZATION_INSTRUMENT
     26#ifndef CLP_FACTORIZATION_NEW_TIMING
     27#define CLP_FACTORIZATION_NEW_TIMING
     28#endif
    2729#ifdef CLP_FACTORIZATION_INSTRUMENT
    2830#include "CoinTime.hpp"
     
    739741          // see if FT
    740742          if (doForrestTomlin_) {
    741                returnCode = CoinFactorization::replaceColumn(regionSparse,
     743             returnCode = CoinFactorization::replaceColumn(regionSparse,
    742744                            pivotRow,
    743745                            pivotCheck,
     
    14041406     }
    14051407}
     1408#ifdef CLP_FACTORIZATION_NEW_TIMING
     1409#ifdef CLP_FACTORIZATION_INSTRUMENT
     1410extern double externalTimeStart;
     1411extern double timeInFactorize;
     1412extern double timeInUpdate;
     1413extern double timeInFactorizeFake;
     1414extern double timeInUpdateFake1;
     1415extern double timeInUpdateFake2;
     1416extern double timeInUpdateTranspose;
     1417extern double timeInUpdateFT;
     1418extern double timeInUpdateTwoFT;
     1419extern double timeInReplace;
     1420extern int numberUpdate;
     1421extern int numberUpdateTranspose;
     1422extern int numberUpdateFT;
     1423extern int numberUpdateTwoFT;
     1424extern int numberReplace;
     1425extern int currentLengthR;
     1426extern int currentLengthU;
     1427extern int currentTakeoutU;
     1428 extern double averageLengthR;
     1429 extern double averageLengthL;
     1430 extern double averageLengthU;
     1431 extern double scaledLengthDense;
     1432 extern double scaledLengthDenseSquared;
     1433 extern double scaledLengthL;
     1434 extern double scaledLengthR;
     1435 extern double scaledLengthU;
     1436extern int startLengthU;
     1437extern int endLengthU;
     1438extern int endLengthU2;
     1439extern int numberAdded;
     1440static int average[3];
     1441static double shortest;
     1442static int ifPrint;
     1443#else
     1444#ifdef CLP_STATIC
     1445static int endLengthU_;
     1446#endif
     1447#endif
     1448#ifdef CLP_STATIC
     1449static double shortestAverage_;
     1450static double totalInR_=0.0;
     1451static double totalInIncreasingU_=0.0;
     1452//static int lastR=0;
     1453//static int lastU=0;
     1454static int lastNumberPivots_=0;
     1455static int effectiveStartNumberU_=0;
     1456#else
     1457//#define shortestAverage shortestAverage_
     1458//#define totalInR totalInR_
     1459//#define totalInIncreasingU totalInIncreasingU_
     1460//#define lastNumberPivots lastNumberPivots_
     1461//#define effectiveStartNumberU effectiveStartNumberU_
     1462//#define endLengthU endLengthU_
     1463#endif
     1464#ifdef CLP_USEFUL_PRINTOUT
     1465static bool readTwiddle=false;
     1466static double weightIncU=1.0;
     1467static double weightR=2.0;
     1468static double weightRest=1.0;
     1469static double weightFactL=30.0;
     1470static double weightFactDense=0.1;
     1471static double weightNrows=10.0;
     1472static double increaseNeeded=1.1;
     1473static double constWeightIterate = 1.0;
     1474static double weightNrowsIterate = 3.0;
     1475#else
     1476#define weightIncU 1.0
     1477#define weightR 2.0
     1478#define weightRest 1.0
     1479#define weightFactL 30.0
     1480#define weightFactDense 0.1
     1481#define weightNrows 10.0
     1482#define increaseNeeded 1.1
     1483#define constWeightIterate   1.0
     1484#define weightNrowsIterate   3.0
     1485#endif
     1486bool
     1487ClpFactorization::timeToRefactorize() const
     1488{
     1489  if (coinFactorizationA_) {
     1490    bool reFactor = (coinFactorizationA_->pivots() * 3 > coinFactorizationA_->maximumPivots() * 2 &&
     1491                     coinFactorizationA_->numberElementsR() * 3 > (coinFactorizationA_->numberElementsL() +
     1492                                                                   coinFactorizationA_->numberElementsU()) * 2 + 1000 &&
     1493                     !coinFactorizationA_->numberDense());
     1494    reFactor=false;
     1495    bool reFactor3=false;
     1496    int numberPivots=coinFactorizationA_->pivots();
     1497    //if (coinFactorizationA_->pivots()<2)
     1498    if (numberPivots>lastNumberPivots_) {
     1499      if (!lastNumberPivots_) {
     1500        //lastR=0;
     1501        //lastU=endLengthU;
     1502        totalInR_=0.0;
     1503        totalInIncreasingU_=0.0;
     1504        shortestAverage_=COIN_DBL_MAX;
     1505#ifdef CLP_USEFUL_PRINTOUT
     1506        if (!readTwiddle) {
     1507          readTwiddle=true;
     1508          char * environ = getenv("CLP_TWIDDLE");
     1509          if (environ) {
     1510            sscanf(environ,"%lg %lg %lg %lg %lg %lg %lg %lg %lg",
     1511      &weightIncU,&weightR,&weightRest,&weightFactL,
     1512      &weightFactDense,&weightNrows,&increaseNeeded,
     1513      &constWeightIterate,&weightNrowsIterate);
     1514          }
     1515            printf("weightIncU %g, weightR %g, weightRest %g, weightFactL %g, weightFactDense %g, weightNrows %g increaseNeeded %g constWeightIterate %g weightNrowsIterate %g\n",
     1516      weightIncU,weightR,weightRest,weightFactL,
     1517        weightFactDense,weightNrows,increaseNeeded,
     1518      constWeightIterate,weightNrowsIterate);
     1519        }
     1520#endif
     1521      }
     1522      lastNumberPivots_=numberPivots;
     1523      int numberDense=coinFactorizationA_->numberDense();
     1524      double nnd=numberDense*numberDense;
     1525      int lengthL=coinFactorizationA_->numberElementsL();
     1526      int lengthR=coinFactorizationA_->numberElementsR();
     1527      int numberRows = coinFactorizationA_->numberRows();
     1528      int lengthU=coinFactorizationA_->numberElementsU()-
     1529        (numberRows-numberDense);
     1530      totalInR_ += lengthR;
     1531      int effectiveU=lengthU-effectiveStartNumberU_;
     1532      totalInIncreasingU_ += effectiveU;
     1533      //lastR=lengthR;
     1534      //lastU=lengthU;
     1535      double rest=lengthL+0.05*nnd;
     1536      double constWeightFactor = weightFactL*lengthL+weightFactDense*nnd
     1537        + weightNrows*numberRows;
     1538      double constWeightIterateX = constWeightIterate*(lengthL+endLengthU_)
     1539        + weightNrowsIterate*numberRows;
     1540      double variableWeight = weightIncU*totalInIncreasingU_+
     1541                               weightR*totalInR_+weightRest*rest;
     1542      double average=constWeightIterateX+
     1543        (constWeightFactor+variableWeight)/static_cast<double>(numberPivots);
     1544#if 0
     1545      if ((numberPivots%20)==0&&!ifPrint3)
     1546      printf("PIV %d nrow %d startU %d now %d L %d R %d dense %g average %g\n",
     1547      numberPivots,numberRows,effectiveStartNumberU_,
     1548      lengthU,lengthL,lengthR,nnd,average);
     1549#endif
     1550      shortestAverage_=CoinMin(shortestAverage_,average);
     1551      if (average>increaseNeeded*shortestAverage_&&
     1552        coinFactorizationA_->pivots()>30) {
     1553        //printf("PIVX %d nrow %d startU %d now %d L %d R %d dense %g average %g\n",
     1554        //numberPivots,numberRows,effectiveStartNumberU_,
     1555        //lengthU,lengthL,lengthR,nnd,average);
     1556        reFactor3=true;
     1557      }
     1558  }
     1559    if (reFactor|| reFactor3) {
     1560#if 0
     1561        printf("%c%cftranCountInput_ %g  ,ftranCountAfterL_ %g  ,ftranCountAfterR_ %g  ,ftranCountAfterU_ %g  ,btranCountInput_ %g  ,btranCountAfterU_ %g  ,btranCountAfterR_ %g  ,btranCountAfterL_ %g  ,numberFtranCounts_ %d  ,numberBtranCounts_ %d  ,ftranAverageAfterL_ %g  ,ftranAverageAfterR_ %g  ,ftranAverageAfterU_ %g  ,btranAverageAfterU_ %g  ,btranAverageAfterR_ %g  ,btranAverageAfterL_ %g\n"
     1562               ,reFactor3 ? 'Y' : 'N'
     1563               ,reFactor ? 'Y' : 'N'
     1564  ,coinFactorizationA_->ftranCountInput_
     1565  ,coinFactorizationA_->ftranCountAfterL_
     1566  ,coinFactorizationA_->ftranCountAfterR_
     1567  ,coinFactorizationA_->ftranCountAfterU_
     1568  ,coinFactorizationA_->btranCountInput_
     1569  ,coinFactorizationA_->btranCountAfterU_
     1570  ,coinFactorizationA_->btranCountAfterR_
     1571  ,coinFactorizationA_->btranCountAfterL_
     1572  ,coinFactorizationA_->numberFtranCounts_
     1573  ,coinFactorizationA_->numberBtranCounts_
     1574  ,coinFactorizationA_->ftranAverageAfterL_
     1575  ,coinFactorizationA_->ftranAverageAfterR_
     1576  ,coinFactorizationA_->ftranAverageAfterU_
     1577  ,coinFactorizationA_->btranAverageAfterU_
     1578  ,coinFactorizationA_->btranAverageAfterR_
     1579               ,coinFactorizationA_->btranAverageAfterL_);
     1580#endif
     1581      reFactor=true;
     1582    }
     1583    return reFactor;
     1584  } else {
     1585    return coinFactorizationB_->pivots() > coinFactorizationB_->numberRows() / 2.45 + 20;
     1586  }
     1587}
     1588#if CLP_FACTORIZATION_NEW_TIMING>1
     1589void
     1590ClpFactorization::statsRefactor(char when) const
     1591{
     1592  int numberPivots=coinFactorizationA_->pivots();
     1593  int numberDense=coinFactorizationA_->numberDense();
     1594  double nnd=numberDense*numberDense;
     1595  int lengthL=coinFactorizationA_->numberElementsL();
     1596  int lengthR=coinFactorizationA_->numberElementsR();
     1597  int numberRows = coinFactorizationA_->numberRows();
     1598  int lengthU=coinFactorizationA_->numberElementsU()-
     1599    (numberRows-numberDense);
     1600  double rest=lengthL+0.05*nnd;
     1601  double constWeightFactor = weightFactL*lengthL+weightFactDense*nnd
     1602    + weightNrows*numberRows;
     1603  double constWeightIterateX = constWeightIterate*(lengthL+endLengthU_)
     1604    + weightNrowsIterate*numberRows;
     1605  double variableWeight = weightIncU*totalInIncreasingU_+
     1606    weightR*totalInR_+weightRest*rest;
     1607  double average=constWeightIterateX+
     1608    (constWeightFactor+variableWeight)/static_cast<double>(numberPivots);
     1609  printf("PIV%c %d nrow %d startU %d now %d L %d R %d dense %g average %g - shortest %g\n",
     1610         when,numberPivots,numberRows,effectiveStartNumberU_,
     1611         lengthU,lengthL,lengthR,nnd,average,shortestAverage_);
     1612}
     1613#endif
     1614#else
     1615bool
     1616ClpFactorization::timeToRefactorize() const
     1617{
     1618    if (coinFactorizationA_) {
     1619    return (coinFactorizationA_->pivots() * 3 > coinFactorizationA_->maximumPivots() * 2 &&
     1620      coinFactorizationA_->numberElementsR() * 3 > (coinFactorizationA_->numberElementsL() +
     1621      coinFactorizationA_->numberElementsU()) * 2 + 1000 &&
     1622      !coinFactorizationA_->numberDense());
     1623  } else {
     1624    return coinFactorizationB_->pivots() > coinFactorizationB_->numberRows() / 2.45 + 20;
     1625  }
     1626#endif
    14061627int
    14071628ClpFactorization::factorize ( ClpSimplex * model,
     
    14201641#ifdef CLP_FACTORIZATION_INSTRUMENT
    14211642     factorization_instrument(-1);
     1643     if (!timeInUpdate) {
     1644      printf("ZZZZ,timeInFactor,nRows,nrows_2,\
     1645startU,endU,lengthL,dense,dense_2,nslacks\n");
     1646      printf("YYYY,time,added,dense,dense_2,averageR,averageL,averageU,\
     1647average2R,average2L,average2U,\
     1648scaledDense,scaledDense_2,scaledL,scaledR,scaledU\n");
     1649      printf("WWWW,time,size0,ratio1,ratio05,ratio02\n");
     1650    }
     1651#endif
     1652#ifdef CLP_FACTORIZATION_INSTRUMENT
     1653     externalTimeStart=CoinCpuTime();
     1654     memset(average,0,sizeof(average));
     1655     shortest=COIN_DBL_MAX;
     1656     ifPrint=0;
     1657     if (!model->numberIterations())
     1658       printf("maxfactor %d\n",coinFactorizationA_->maximumPivots());
    14221659#endif
    14231660     bool anyChanged = false;
     
    19992236                    else
    20002237                         numberElements = 0;
     2238#ifdef CLP_FACTORIZATION_NEW_TIMING
     2239                    lastNumberPivots_=0;
     2240                    effectiveStartNumberU_=numberElements-numberRows;
     2241                    //printf("%d slacks,%d in U at beginning\n",
     2242                    //numberRowBasic,numberElements);
     2243#endif
    20012244                    coinFactorizationA_->setNumberElementsU(numberElements);
    20022245                    //saveFactorization("dump.d");
     
    20062249                         coinFactorizationA_->preProcess ( 3 ); // no row copy
    20072250                    coinFactorizationA_->factor (  );
     2251#ifdef CLP_FACTORIZATION_NEW_TIMING
     2252                    endLengthU_ = coinFactorizationA_->numberElements() -
     2253                      coinFactorizationA_->numberDense()*coinFactorizationA_->numberDense()
     2254                      -coinFactorizationA_->numberElementsL();
     2255#endif
    20082256                    if (coinFactorizationA_->status() == -99) {
    20092257                         // get more memory
     
    21082356               } else if (coinFactorizationA_->status() == -1 && (solveType == 0 || solveType == 2)) {
    21092357                    // This needs redoing as it was merged coding - does not need array
     2358#if 1
    21102359                    int numberTotal = numberRows + numberColumns;
    21112360                    int * isBasic = new int [numberTotal];
     
    21682417                    }
    21692418                    delete [] isBasic;
     2419#else
     2420                    {
     2421                      //int * lastColumn = lastColumn_.array(); // -1 or pivot row
     2422                      int * lastRow = coinFactorizationA_->lastRow(); // -1 or pivot sequence (inside sequence)
     2423                      for ( int i=0;i<numberRows;i++) {
     2424                        int iSeq = lastRow[i];
     2425                        if (iSeq >=0) {
     2426                          pivotVariable[i]=pivotTemp[iSeq];
     2427                          model->setRowStatus(i, ClpSimplex::superBasic);
     2428                        } else {
     2429                          pivotVariable[i]=i+numberColumns;
     2430                        }
     2431                      }
     2432                    }
     2433#endif
    21702434                    double * columnLower = model->lowerRegion();
    21712435                    double * columnUpper = model->upperRegion();
     
    23042568     if (!networkBasis_) {
    23052569#endif
     2570#ifdef CLP_FACTORIZATION_NEW_TIMING
    23062571#ifdef CLP_FACTORIZATION_INSTRUMENT
    23072572          factorization_instrument(-1);
     2573#endif
     2574          int nOld=0;
     2575          int nNew=0;
     2576          int seq;
     2577          const CoinPackedMatrix * matrix=model->matrix();
     2578          const int * columnLength = matrix->getVectorLengths();
     2579          seq=model->sequenceIn();
     2580          if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
     2581            if (seq<model->numberColumns())
     2582              nNew=columnLength[seq];
     2583            else
     2584              nNew=1;
     2585          }
     2586          seq=model->sequenceOut();
     2587          if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
     2588            if (seq<model->numberColumns())
     2589              nOld=columnLength[seq];
     2590            else
     2591              nOld=1;
     2592          }
     2593          effectiveStartNumberU_ += nNew-nOld;
    23082594#endif
    23092595          int returnCode;
     
    23112597          if (!coinFactorizationA_ || coinFactorizationA_->forrestTomlin()) {
    23122598               if (coinFactorizationA_) {
    2313                     returnCode = coinFactorizationA_->replaceColumn(regionSparse,
    2314                                  pivotRow,
    2315                                  pivotCheck,
    2316                                  checkBeforeModifying,
    2317                                  acceptablePivot);
     2599#if ABC_USE_COIN_FACTORIZATION<2
     2600                 returnCode =
     2601                   coinFactorizationA_->replaceColumn(regionSparse,
     2602                                                      pivotRow,
     2603                                                      pivotCheck,
     2604                                                      checkBeforeModifying,
     2605                                                      acceptablePivot);
     2606#else
     2607                 // fake btran alpha until I understand
     2608                 double btranAlpha=model->alpha();
     2609                 double ftAlpha =
     2610                   coinFactorizationA_->checkReplacePart1(regionSparse,
     2611                                                          pivotRow);
     2612                 returnCode =
     2613                   coinFactorizationA_->checkReplacePart2(pivotRow,
     2614                                                          btranAlpha,
     2615                                                          model->alpha(),
     2616                                                          ftAlpha,
     2617                                                          acceptablePivot);
     2618                 if (returnCode<2)
     2619                   coinFactorizationA_->replaceColumnPart3(regionSparse,
     2620                                                           pivotRow,
     2621                                                           model->alpha());
     2622#endif
    23182623               } else {
    23192624                    bool tab = coinFactorizationB_->wantsTableauColumn();
  • trunk/Clp/src/ClpFactorization.hpp

    r1665 r2078  
    2424#ifndef COIN_FAST_CODE
    2525#define COIN_FAST_CODE
     26#endif
     27#ifndef CLP_FACTORIZATION_NEW_TIMING
     28#define CLP_FACTORIZATION_NEW_TIMING 1
    2629#endif
    2730
     
    222225     }
    223226#endif
    224      inline bool timeToRefactorize() const {
    225           if (coinFactorizationA_) {
    226                return (coinFactorizationA_->pivots() * 3 > coinFactorizationA_->maximumPivots() * 2 &&
    227                        coinFactorizationA_->numberElementsR() * 3 > (coinFactorizationA_->numberElementsL() +
    228                                  coinFactorizationA_->numberElementsU()) * 2 + 1000 &&
    229                        !coinFactorizationA_->numberDense());
    230           } else {
    231                return coinFactorizationB_->pivots() > coinFactorizationB_->numberRows() / 2.45 + 20;
    232           }
    233      }
     227     bool timeToRefactorize() const;
     228#if CLP_FACTORIZATION_NEW_TIMING>1
     229     void statsRefactor(char when) const;
     230#endif
    234231     /// Level of detail of messages
    235232     inline int messageLevel (  ) const {
     
    421418     int goDenseThreshold_;
    422419#endif
     420#ifdef CLP_FACTORIZATION_NEW_TIMING
     421     /// For guessing when to re-factorize
     422     mutable double shortestAverage_;
     423     mutable double totalInR_;
     424     mutable double totalInIncreasingU_;
     425     mutable int endLengthU_;
     426     mutable int lastNumberPivots_;
     427     mutable int effectiveStartNumberU_;
     428#endif
    423429     //@}
    424430};
    425 
    426 #endif
     431 
     432#endif
  • trunk/Clp/src/ClpMain.cpp

    r2074 r2078  
    169169#endif
    170170#endif
     171//#define CLP_MALLOC_STATISTICS
     172#ifdef CLP_MALLOC_STATISTICS
     173#include <malloc.h>
     174#include <exception>
     175#include <new>
     176static double malloc_times = 0.0;
     177static double malloc_total = 0.0;
     178static int malloc_amount[] = {0, 32, 128, 256, 1024, 4096, 16384, 65536, 262144, COIN_INT_MAX};
     179static int malloc_n = 10;
     180double malloc_counts[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
     181bool malloc_counts_on = true;
     182void * operator new (size_t size) throw (std::bad_alloc)
     183{
     184    malloc_times ++;
     185    malloc_total += size;
     186    int i;
     187    for (i = 0; i < malloc_n; i++) {
     188        if ((int) size <= malloc_amount[i]) {
     189            malloc_counts[i]++;
     190            break;
     191        }
     192    }
     193    if (size>100000) {
     194      printf("allocating %ld bytes\n",size);
     195    }
     196    void * p = malloc(size);
     197    return p;
     198}
     199void operator delete (void *p) throw()
     200{
     201    free(p);
     202}
     203static void malloc_stats2()
     204{
     205    double average = malloc_total / malloc_times;
     206    printf("count %g bytes %g - average %g\n", malloc_times, malloc_total, average);
     207    for (int i = 0; i < malloc_n; i++)
     208        printf("%g ", malloc_counts[i]);
     209    printf("\n");
     210    malloc_times = 0.0;
     211    malloc_total = 0.0;
     212    memset(malloc_counts, 0, sizeof(malloc_counts));
     213    // print results
     214}
     215#endif  //CLP_MALLOC_STATISTICS
    171216int
    172217#if defined(_MSC_VER)
  • trunk/Clp/src/ClpMatrixBase.hpp

    r1722 r2078  
    1414class ClpSimplex;
    1515class ClpModel;
     16// Compilers can produce better code if they know about __restrict
     17#ifndef COIN_RESTRICT
     18#ifdef COIN_USE_RESTRICT
     19#define COIN_RESTRICT __restrict
     20#else
     21#define COIN_RESTRICT
     22#endif
     23#endif
    1624
    1725/** Abstract base class for Clp Matrices
     
    264272         @pre <code>y</code> must be of size <code>numRows()</code> */
    265273     virtual void times(double scalar,
    266                         const double * x, double * y) const = 0;
     274                        const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const = 0;
    267275     /** And for scaling - default aborts for when scaling not supported
    268276         (unless pointers NULL when as normal)
    269277     */
    270278     virtual void times(double scalar,
    271                         const double * x, double * y,
    272                         const double * rowScale,
    273                         const double * columnScale) const;
     279                        const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
     280                        const double * COIN_RESTRICT rowScale,
     281                        const double * COIN_RESTRICT columnScale) const;
    274282     /** Return <code>y + x * scalar * A</code> in <code>y</code>.
    275283         @pre <code>x</code> must be of size <code>numRows()</code>
    276284         @pre <code>y</code> must be of size <code>numColumns()</code> */
    277285     virtual void transposeTimes(double scalar,
    278                                  const double * x, double * y) const = 0;
     286                                 const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const = 0;
    279287     /** And for scaling - default aborts for when scaling not supported
    280288         (unless pointers NULL when as normal)
    281289     */
    282290     virtual void transposeTimes(double scalar,
    283                                  const double * x, double * y,
    284                                  const double * rowScale,
    285                                  const double * columnScale,
    286                                  double * spare = NULL) const;
     291                                 const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
     292                                 const double * COIN_RESTRICT rowScale,
     293                                 const double * COIN_RESTRICT columnScale,
     294                                 double * COIN_RESTRICT spare = NULL) const;
    287295#if COIN_LONG_WORK
    288296     // For long double versions (aborts if not supported)
    289297     virtual void times(CoinWorkDouble scalar,
    290                         const CoinWorkDouble * x, CoinWorkDouble * y) const ;
     298                        const CoinWorkDouble * COIN_RESTRICT x, CoinWorkDouble * COIN_RESTRICT y) const ;
    291299     virtual void transposeTimes(CoinWorkDouble scalar,
    292                                  const CoinWorkDouble * x, CoinWorkDouble * y) const ;
     300                                 const CoinWorkDouble * COIN_RESTRICT x, CoinWorkDouble * COIN_RESTRICT y) const ;
    293301#endif
    294302     /** Return <code>x * scalar *A + y</code> in <code>z</code>.
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r2074 r2078  
    276276void
    277277ClpPackedMatrix::times(double scalar,
    278                        const double * x, double * y) const
     278                       const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
    279279{
    280280     int iRow, iColumn;
    281281     // get matrix data pointers
    282      const int * row = matrix_->getIndices();
    283      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    284      const double * elementByColumn = matrix_->getElements();
     282     const int * COIN_RESTRICT row = matrix_->getIndices();
     283     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     284     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    285285     //memset(y,0,matrix_->getNumRows()*sizeof(double));
    286286     assert (((flags_&0x02) != 0) == matrix_->hasGaps());
     
    318318void
    319319ClpPackedMatrix::transposeTimes(double scalar,
    320                                 const double * x, double * y) const
     320                                const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
    321321{
    322322     int iColumn;
    323323     // get matrix data pointers
    324      const int * row = matrix_->getIndices();
    325      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    326      const double * elementByColumn = matrix_->getElements();
     324     const int * COIN_RESTRICT row = matrix_->getIndices();
     325     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     326     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    327327     if (!(flags_ & 2)) {
    328328          if (scalar == -1.0) {
     
    593593{
    594594     columnArray->clear();
    595      double * pi = rowArray->denseVector();
     595     double * COIN_RESTRICT pi = rowArray->denseVector();
    596596     int numberNonZero = 0;
    597      int * index = columnArray->getIndices();
    598      double * array = columnArray->denseVector();
     597     int * COIN_RESTRICT index = columnArray->getIndices();
     598     double * COIN_RESTRICT array = columnArray->denseVector();
    599599     int numberInRowArray = rowArray->getNumElements();
    600600     // maybe I need one in OsiSimplex
     
    645645          int iColumn;
    646646          // get matrix data pointers
    647           const int * row = matrix_->getIndices();
    648           const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    649           const int * columnLength = matrix_->getVectorLengths();
    650           const double * elementByColumn = matrix_->getElements();
    651           const double * rowScale = model->rowScale();
     647          const int * COIN_RESTRICT row = matrix_->getIndices();
     648          const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     649          const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     650          const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     651          const double * COIN_RESTRICT rowScale = model->rowScale();
    652652#if 0
    653653          ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
     
    664664               // need to expand pi into y
    665665               assert(y->capacity() >= numberRows);
    666                double * piOld = pi;
     666               double * COIN_RESTRICT piOld = pi;
    667667               pi = y->denseVector();
    668                const int * whichRow = rowArray->getIndices();
     668               const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    669669               int i;
    670670               if (!rowScale) {
     
    884884          // need to expand pi into y
    885885          assert(y->capacity() >= model->numberRows());
    886           double * piOld = pi;
     886          double * COIN_RESTRICT piOld = pi;
    887887          pi = y->denseVector();
    888888          const int * COIN_RESTRICT whichRow = rowArray->getIndices();
     
    907907                              CoinIndexedVector * spareArray = model->rowArray(3);
    908908                              // also do dualColumn stuff
    909                               double * spare = spareArray->denseVector();
    910                               int * spareIndex = spareArray->getIndices();
    911                               const double * reducedCost = model->djRegion(0);
     909                              double * COIN_RESTRICT spare = spareArray->denseVector();
     910                              int * COIN_RESTRICT spareIndex = spareArray->getIndices();
     911                              const double * COIN_RESTRICT reducedCost = model->djRegion(0);
    912912                              double multiplier[] = { -1.0, 1.0};
    913913                              double dualT = - model->currentDualTolerance();
     
    918918                              double bestPossible = 0.0;
    919919                              int addSequence = model->numberColumns();
    920                               const unsigned char * statusArray = model->statusArray() + addSequence;
     920                              const unsigned char * COIN_RESTRICT statusArray = model->statusArray() + addSequence;
    921921                              int numberRemaining = 0;
    922922                              assert (scalar == -1.0);
     
    11061106               // scaled
    11071107               if (scalar == -1.0) {
    1108                     const double * columnScale = model->columnScale();
     1108                    const double * COIN_RESTRICT columnScale = model->columnScale();
    11091109                    double value = 0.0;
    11101110                    double scale = columnScale[0];
     
    11361136                    }
    11371137               } else {
    1138                     const double * columnScale = model->columnScale();
     1138                    const double * COIN_RESTRICT columnScale = model->columnScale();
    11391139                    double value = 0.0;
    11401140                    double scale = columnScale[0] * scalar;
     
    11821182{
    11831183     columnArray->clear();
    1184      double * pi = rowArray->denseVector();
     1184     double * COIN_RESTRICT pi = rowArray->denseVector();
    11851185     int numberNonZero = 0;
    1186      int * index = columnArray->getIndices();
    1187      double * array = columnArray->denseVector();
     1186     int * COIN_RESTRICT index = columnArray->getIndices();
     1187     double * COIN_RESTRICT array = columnArray->denseVector();
    11881188     int numberInRowArray = rowArray->getNumElements();
    11891189     // maybe I need one in OsiSimplex
    11901190     double zeroTolerance = model->zeroTolerance();
    1191      const int * column = matrix_->getIndices();
    1192      const CoinBigIndex * rowStart = getVectorStarts();
    1193      const double * element = getElements();
    1194      const int * whichRow = rowArray->getIndices();
     1191     const int * COIN_RESTRICT column = matrix_->getIndices();
     1192     const CoinBigIndex * COIN_RESTRICT rowStart = getVectorStarts();
     1193     const double * COIN_RESTRICT element = getElements();
     1194     const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    11951195     bool packed = rowArray->packedMode();
    11961196     if (numberInRowArray > 2) {
     
    12011201          int numberOriginal = 0;
    12021202          if (packed) {
    1203                int * index = columnArray->getIndices();
    1204                double * array = columnArray->denseVector();
     1203               int * COIN_RESTRICT index = columnArray->getIndices();
     1204               double * COIN_RESTRICT array = columnArray->denseVector();
    12051205#if 0
    12061206               {
    1207                  double  * array2 = y->denseVector();
     1207                 double  * COIN_RESTRICT array2 = y->denseVector();
    12081208                 int numberColumns = matrix_->getNumCols();
    12091209                 for (int i=0;i<numberColumns;i++) {
     
    12301230#if COIN_SPARSE_MATRIX != 2
    12311231                 // and set up mark as char array
    1232                  char * marked = reinterpret_cast<char *> (index+columnArray->capacity());
    1233                  int * lookup = y->getIndices();
     1232                 char * COIN_RESTRICT marked = reinterpret_cast<char *> (index+columnArray->capacity());
     1233                 int * COIN_RESTRICT lookup = y->getIndices();
    12341234#ifndef NDEBUG
    12351235                 //int numberColumns = matrix_->getNumCols();
     
    12401240                                                             lookup,marked,zeroTolerance,scalar);
    12411241#else
    1242                  double  * array2 = y->denseVector();
     1242                 double  * COIN_RESTRICT array2 = y->denseVector();
    12431243                 numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
    12441244                                                            array2,zeroTolerance,scalar);
     
    12501250               columnArray->setNumElements(numberNonZero);
    12511251          } else {
    1252                double * markVector = y->denseVector();
     1252               double * COIN_RESTRICT markVector = y->denseVector();
    12531253               numberNonZero = 0;
    12541254               // and set up mark as char array
    1255                char * marked = reinterpret_cast<char *> (markVector);
     1255               char * COIN_RESTRICT marked = reinterpret_cast<char *> (markVector);
    12561256               for (i = 0; i < numberOriginal; i++) {
    12571257                    int iColumn = index[i];
     
    17621762          CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
    17631763{
    1764      double * pi = piVector->denseVector();
     1764     double * COIN_RESTRICT pi = piVector->denseVector();
    17651765     int numberNonZero = 0;
    1766      int * index = output->getIndices();
    1767      double * array = output->denseVector();
    1768      const int * column = matrix_->getIndices();
    1769      const CoinBigIndex * rowStart = matrix_->getVectorStarts();
    1770      const double * element = matrix_->getElements();
    1771      const int * whichRow = piVector->getIndices();
     1766     int * COIN_RESTRICT index = output->getIndices();
     1767     double * COIN_RESTRICT array = output->denseVector();
     1768     const int * COIN_RESTRICT column = matrix_->getIndices();
     1769     const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     1770     const double * COIN_RESTRICT element = matrix_->getElements();
     1771     const int * COIN_RESTRICT whichRow = piVector->getIndices();
    17721772     int iRow0 = whichRow[0];
    17731773     int iRow1 = whichRow[1];
     
    17831783     }
    17841784     // and set up mark as char array
    1785      char * marked = reinterpret_cast<char *> (index + output->capacity());
    1786      int * lookup = spareVector->getIndices();
     1785     char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + output->capacity());
     1786     int * COIN_RESTRICT lookup = spareVector->getIndices();
    17871787     double value = pi0 * scalar;
    17881788     CoinBigIndex j;
     
    18601860          const double tolerance, const double scalar) const
    18611861{
    1862      double * pi = piVector->denseVector();
     1862     double * COIN_RESTRICT pi = piVector->denseVector();
    18631863     int numberNonZero = 0;
    1864      int * index = output->getIndices();
    1865      double * array = output->denseVector();
    1866      const int * column = matrix_->getIndices();
    1867      const CoinBigIndex * rowStart = matrix_->getVectorStarts();
    1868      const double * element = matrix_->getElements();
     1864     int * COIN_RESTRICT index = output->getIndices();
     1865     double * COIN_RESTRICT array = output->denseVector();
     1866     const int * COIN_RESTRICT column = matrix_->getIndices();
     1867     const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     1868     const double * COIN_RESTRICT element = matrix_->getElements();
    18691869     int iRow = piVector->getIndices()[0];
    18701870     numberNonZero = 0;
     
    20632063                                 double referenceIn, double devex,
    20642064                                 // Array for exact devex to say what is in reference framework
    2065                                  unsigned int * reference,
    2066                                  double * weights, double scaleFactor)
     2065                                 unsigned int * COIN_RESTRICT reference,
     2066                                 double * COIN_RESTRICT weights, double scaleFactor)
    20672067{
    20682068     // put row of tableau in dj1
    2069      double * pi = pi1->denseVector();
     2069     double * COIN_RESTRICT pi = pi1->denseVector();
    20702070     int numberNonZero = 0;
    2071      int * index = dj1->getIndices();
    2072      double * array = dj1->denseVector();
     2071     int * COIN_RESTRICT index = dj1->getIndices();
     2072     double * COIN_RESTRICT array = dj1->denseVector();
    20732073     int numberInRowArray = pi1->getNumElements();
    20742074     double zeroTolerance = model->zeroTolerance();
     
    20772077     int iColumn;
    20782078     // get matrix data pointers
    2079      const int * row = matrix_->getIndices();
    2080      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    2081      const double * elementByColumn = matrix_->getElements();
    2082      const double * rowScale = model->rowScale();
     2079     const int * COIN_RESTRICT row = matrix_->getIndices();
     2080     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     2081     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     2082     const double * COIN_RESTRICT rowScale = model->rowScale();
    20832083     assert (!spare->getNumElements());
    20842084     assert (numberActiveColumns_ > 0);
    2085      double * piWeight = pi2->denseVector();
     2085     double * COIN_RESTRICT piWeight = pi2->denseVector();
    20862086     assert (!pi2->packedMode());
    20872087     bool killDjs = (scaleFactor == 0.0);
     
    20912091          // need to expand pi into y
    20922092          assert(spare->capacity() >= model->numberRows());
    2093           double * piOld = pi;
     2093          double * COIN_RESTRICT piOld = pi;
    20942094          pi = spare->denseVector();
    2095           const int * whichRow = pi1->getIndices();
     2095          const int * COIN_RESTRICT whichRow = pi1->getIndices();
    20962096          int i;
    20972097          ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
     
    21762176               }
    21772177               if (!columnCopy_) {
    2178                     const double * columnScale = model->columnScale();
     2178                    const double * COIN_RESTRICT columnScale = model->columnScale();
    21792179                    CoinBigIndex j;
    21802180                    CoinBigIndex end = columnStart[0];
     
    22932293               // can also scale piWeight as not used again
    22942294               int numberWeight = pi2->getNumElements();
    2295                const int * indexWeight = pi2->getIndices();
     2295               const int * COIN_RESTRICT indexWeight = pi2->getIndices();
    22962296               for (int i = 0; i < numberWeight; i++) {
    22972297                    int iRow = indexWeight[i];
    22982298                    piWeight[iRow] *= rowScale[iRow];
    22992299               }
    2300                const double * columnScale = model->columnScale();
     2300               const double * COIN_RESTRICT columnScale = model->columnScale();
    23012301               CoinBigIndex j;
    23022302               CoinBigIndex end = columnStart[0];
     
    23572357                              double referenceIn, double devex,
    23582358                              // Array for exact devex to say what is in reference framework
    2359                               unsigned int * reference,
    2360                               double * weights, double scaleFactor)
     2359                              unsigned int * COIN_RESTRICT reference,
     2360                              double * COIN_RESTRICT weights, double scaleFactor)
    23612361{
    23622362     int number = dj1->getNumElements();
    2363      const int * index = dj1->getIndices();
    2364      double * array = dj1->denseVector();
     2363     const int * COIN_RESTRICT index = dj1->getIndices();
     2364     double * COIN_RESTRICT array = dj1->denseVector();
    23652365     assert( dj1->packedMode());
    23662366
    23672367     // get matrix data pointers
    2368      const int * row = matrix_->getIndices();
    2369      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    2370      const int * columnLength = matrix_->getVectorLengths();
    2371      const double * elementByColumn = matrix_->getElements();
    2372      const double * rowScale = model->rowScale();
    2373      double * piWeight = pi2->denseVector();
     2368     const int * COIN_RESTRICT row = matrix_->getIndices();
     2369     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     2370     const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     2371     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     2372     const double * COIN_RESTRICT rowScale = model->rowScale();
     2373     double * COIN_RESTRICT piWeight = pi2->denseVector();
    23742374     bool killDjs = (scaleFactor == 0.0);
    23752375     if (!scaleFactor)
     
    32883288#endif
    32893289     ClpMatrixBase * rowCopyBase = model->rowCopy();
    3290      double * rowScale;
    3291      double * columnScale;
     3290     double * COIN_RESTRICT rowScale;
     3291     double * COIN_RESTRICT columnScale;
    32923292     //assert (!model->rowScale());
    32933293     bool arraysExist;
    3294      double * inverseRowScale = NULL;
    3295      double * inverseColumnScale = NULL;
     3294     double * COIN_RESTRICT inverseRowScale = NULL;
     3295     double * COIN_RESTRICT inverseColumnScale = NULL;
    32963296     if (!model->rowScale()) {
    32973297          rowScale = new double [numberRows*2];
     
    33103310     assert (inverseColumnScale == columnScale + numberColumns);
    33113311     // we are going to mark bits we are interested in
    3312      char * usefulColumn = new char [numberColumns];
    3313      double * rowLower = model->rowLower();
    3314      double * rowUpper = model->rowUpper();
    3315      double * columnLower = model->columnLower();
    3316      double * columnUpper = model->columnUpper();
     3312     char * COIN_RESTRICT usefulColumn = new char [numberColumns];
     3313     double * COIN_RESTRICT rowLower = model->rowLower();
     3314     double * COIN_RESTRICT rowUpper = model->rowUpper();
     3315     double * COIN_RESTRICT columnLower = model->columnLower();
     3316     double * COIN_RESTRICT columnUpper = model->columnUpper();
    33173317     int iColumn, iRow;
    33183318     //#define LEAVE_FIXED
     
    33243324     double smallest = 1.0e50;
    33253325     // get matrix data pointers
    3326      int * row = matrix_->getMutableIndices();
    3327      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    3328      int * columnLength = matrix_->getMutableVectorLengths();
    3329      double * elementByColumn = matrix_->getMutableElements();
     3326     int * COIN_RESTRICT row = matrix_->getMutableIndices();
     3327     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     3328     int * COIN_RESTRICT columnLength = matrix_->getMutableVectorLengths();
     3329     double * COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
    33303330     int deletedElements = 0;
    33313331     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     
    33793379     // don't scale integers if option set
    33803380     if ((model->specialOptions()&4194304)!=0 && model->integerInformation()) {
    3381        const char * integer  = model->integerInformation();
     3381       const char * COIN_RESTRICT integer  = model->integerInformation();
    33823382       for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    33833383         if (integer[iColumn])
     
    34273427#endif
    34283428
    3429           const int * column = rowCopy->getIndices();
    3430           const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    3431           const double * element = rowCopy->getElements();
     3429          const int * COIN_RESTRICT column = rowCopy->getIndices();
     3430          const CoinBigIndex * COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
     3431          const double * COIN_RESTRICT element = rowCopy->getElements();
    34323432          // need to scale
    34333433          if (largest > 1.0e13 * smallest) {
     
    34933493#ifdef USE_OBJECTIVE
    34943494                    // This will be used to help get scale factors
    3495                     double * objective = new double[numberColumns];
     3495                    double * COIN_RESTRICT objective = new double[numberColumns];
    34963496                    CoinMemcpyN(model->costRegion(1), numberColumns, objective);
    34973497                    double objScale = 1.0;
     
    37703770                    static_cast< ClpPackedMatrix*>(rowCopyBase);
    37713771#endif
    3772                const int * column = rowCopy->getIndices();
    3773                const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     3772               const int * COIN_RESTRICT column = rowCopy->getIndices();
     3773               const CoinBigIndex * COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
    37743774               CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    37753775               int numberXColumns = quadratic->getNumCols();
     
    37803780                    //const int * columnQuadratic = quadratic->getIndices();
    37813781                    //const int * columnQuadraticStart = quadratic->getVectorStarts();
    3782                     const int * columnQuadraticLength = quadratic->getVectorLengths();
     3782                    const int * COIN_RESTRICT columnQuadraticLength = quadratic->getVectorLengths();
    37833783                    for (i = 0; i < numberXColumns; i++) {
    37843784                         int length = columnQuadraticLength[i];
     
    38443844               ClpPackedMatrix* rowCopy =
    38453845                    static_cast< ClpPackedMatrix*>(model->rowCopy());
    3846                double * element = rowCopy->getMutableElements();
    3847                const int * column = rowCopy->getIndices();
    3848                const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     3846               double * COIN_RESTRICT element = rowCopy->getMutableElements();
     3847               const int * COIN_RESTRICT column = rowCopy->getIndices();
     3848               const CoinBigIndex * COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
    38493849               // scale row copy
    38503850               for (iRow = 0; iRow < numberRows; iRow++) {
    38513851                    CoinBigIndex j;
    38523852                    double scale = rowScale[iRow];
    3853                     double * elementsInThisRow = element + rowStart[iRow];
    3854                     const int * columnsInThisRow = column + rowStart[iRow];
     3853                    double * COIN_RESTRICT elementsInThisRow = element + rowStart[iRow];
     3854                    const int * COIN_RESTRICT columnsInThisRow = column + rowStart[iRow];
    38553855                    int number = rowStart[iRow+1] - rowStart[iRow];
    38563856                    assert (number <= numberColumns);
     
    38683868                    model->setClpScaledMatrix(scaled);
    38693869                    // get matrix data pointers
    3870                     const int * row = scaledMatrix->getIndices();
    3871                     const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
     3870                    const int * COIN_RESTRICT row = scaledMatrix->getIndices();
     3871                    const CoinBigIndex * COIN_RESTRICT columnStart = scaledMatrix->getVectorStarts();
    38723872#ifndef NDEBUG
    3873                     const int * columnLength = scaledMatrix->getVectorLengths();
    3874 #endif
    3875                     double * elementByColumn = scaledMatrix->getMutableElements();
     3873                    const int * COIN_RESTRICT columnLength = scaledMatrix->getVectorLengths();
     3874#endif
     3875                    double * COIN_RESTRICT elementByColumn = scaledMatrix->getMutableElements();
    38763876                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    38773877                         CoinBigIndex j;
     
    39093909     if (!model->rowScale())
    39103910          return;
    3911      double * rowScale = model->mutableRowScale();
    3912      double * columnScale = model->mutableColumnScale();
     3911     double * COIN_RESTRICT rowScale = model->mutableRowScale();
     3912     double * COIN_RESTRICT columnScale = model->mutableColumnScale();
    39133913     // copy without gaps
    39143914     CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, 0, 0);
     
    39163916     model->setClpScaledMatrix(scaled);
    39173917     // get matrix data pointers
    3918      const int * row = scaledMatrix->getIndices();
    3919      const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
     3918     const int * COIN_RESTRICT row = scaledMatrix->getIndices();
     3919     const CoinBigIndex * COIN_RESTRICT columnStart = scaledMatrix->getVectorStarts();
    39203920#ifndef NDEBUG
    3921      const int * columnLength = scaledMatrix->getVectorLengths();
    3922 #endif
    3923      double * elementByColumn = scaledMatrix->getMutableElements();
     3921     const int * COIN_RESTRICT columnLength = scaledMatrix->getVectorLengths();
     3922#endif
     3923     double * COIN_RESTRICT elementByColumn = scaledMatrix->getMutableElements();
    39243924     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    39253925          CoinBigIndex j;
     
    39423942                        int iColumn) const
    39433943{
    3944      const double * rowScale = model->rowScale();
    3945      const int * row = matrix_->getIndices();
    3946      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    3947      const int * columnLength = matrix_->getVectorLengths();
    3948      const double * elementByColumn = matrix_->getElements();
     3944     const double * COIN_RESTRICT rowScale = model->rowScale();
     3945     const int * COIN_RESTRICT row = matrix_->getIndices();
     3946     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     3947     const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     3948     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    39493949     CoinBigIndex i;
    39503950     if (!rowScale) {
     
    39723972                              int iColumn) const
    39733973{
    3974      const double * rowScale = model->rowScale();
    3975      const int * row = matrix_->getIndices();
    3976      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    3977      const int * columnLength = matrix_->getVectorLengths();
    3978      const double * elementByColumn = matrix_->getElements();
     3974     const double * COIN_RESTRICT rowScale = model->rowScale();
     3975     const int * COIN_RESTRICT row = matrix_->getIndices();
     3976     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     3977     const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     3978     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    39793979     CoinBigIndex i;
    3980      int * index = rowArray->getIndices();
    3981      double * array = rowArray->denseVector();
     3980     int * COIN_RESTRICT index = rowArray->getIndices();
     3981     double * COIN_RESTRICT array = rowArray->denseVector();
    39823982     int number = 0;
    39833983     if (!rowScale) {
     
    40154015                     int iColumn, double multiplier) const
    40164016{
    4017      const double * rowScale = model->rowScale();
    4018      const int * row = matrix_->getIndices();
    4019      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    4020      const int * columnLength = matrix_->getVectorLengths();
    4021      const double * elementByColumn = matrix_->getElements();
     4017     const double * COIN_RESTRICT rowScale = model->rowScale();
     4018     const int * COIN_RESTRICT row = matrix_->getIndices();
     4019     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     4020     const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     4021     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    40224022     CoinBigIndex i;
    40234023     if (!rowScale) {
     
    40394039/* Adds multiple of a column into an array */
    40404040void
    4041 ClpPackedMatrix::add(const ClpSimplex * model, double * array,
     4041ClpPackedMatrix::add(const ClpSimplex * model, double * COIN_RESTRICT array,
    40424042                     int iColumn, double multiplier) const
    40434043{
    4044      const double * rowScale = model->rowScale();
    4045      const int * row = matrix_->getIndices();
    4046      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    4047      const int * columnLength = matrix_->getVectorLengths();
    4048      const double * elementByColumn = matrix_->getElements();
     4044     const double * COIN_RESTRICT rowScale = model->rowScale();
     4045     const int * COIN_RESTRICT row = matrix_->getIndices();
     4046     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     4047     const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     4048     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    40494049     CoinBigIndex i;
    40504050     if (!rowScale) {
     
    40934093     double firstBadElement = 0.0;
    40944094     // get matrix data pointers
    4095      const int * row = matrix_->getIndices();
    4096      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    4097      const int * columnLength = matrix_->getVectorLengths();
    4098      const double * elementByColumn = matrix_->getElements();
     4095     const int * COIN_RESTRICT row = matrix_->getIndices();
     4096     const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     4097     const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     4098     const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    40994099     int numberRows = model->numberRows();
    41004100     int numberColumns = matrix_->getNumCols();
     
    41254125     assert (check == 15 || check == 11);
    41264126     if (check == 15) {
    4127           int * mark = new int [numberRows];
     4127          int * COIN_RESTRICT mark = new int [numberRows];
    41284128          int i;
    41294129          for (i = 0; i < numberRows; i++)
     
    42584258     const double * COIN_RESTRICT element = matrix_->getElements();
    42594259     const int * COIN_RESTRICT whichRow = piVector->getIndices();
    4260      int * fakeRow = const_cast<int *> (whichRow);
     4260     int * COIN_RESTRICT fakeRow = const_cast<int *> (whichRow);
    42614261     fakeRow[numberInRowArray]=0; // so can touch
    42624262#ifndef NDEBUG
     
    44544454     int start = static_cast<int> (startFraction * numberActiveColumns_);
    44554455     int end = CoinMin(static_cast<int> (endFraction * numberActiveColumns_ + 1), numberActiveColumns_);
    4456      const double * element = matrix_->getElements();
    4457      const int * row = matrix_->getIndices();
    4458      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    4459      const int * length = matrix_->getVectorLengths();
    4460      const double * rowScale = model->rowScale();
    4461      const double * columnScale = model->columnScale();
     4456     const double * COIN_RESTRICT element = matrix_->getElements();
     4457     const int * COIN_RESTRICT row = matrix_->getIndices();
     4458     const CoinBigIndex * COIN_RESTRICT startColumn = matrix_->getVectorStarts();
     4459     const int * COIN_RESTRICT length = matrix_->getVectorLengths();
     4460     const double * COIN_RESTRICT rowScale = model->rowScale();
     4461     const double * COIN_RESTRICT columnScale = model->columnScale();
    44624462     int iSequence;
    44634463     CoinBigIndex j;
    44644464     double tolerance = model->currentDualTolerance();
    4465      double * reducedCost = model->djRegion();
    4466      const double * duals = model->dualRowSolution();
    4467      const double * cost = model->costRegion();
     4465     double * COIN_RESTRICT reducedCost = model->djRegion();
     4466     const double * COIN_RESTRICT duals = model->dualRowSolution();
     4467     const double * COIN_RESTRICT cost = model->costRegion();
    44684468     double bestDj;
    44694469     if (bestSequence >= 0)
     
    47614761#endif
    47624762
    4763           const int * column = rowCopy->getIndices();
    4764           const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    4765           double * element = rowCopy->getMutableElements();
    4766           const double * rowScale = model->rowScale();
    4767           const double * columnScale = model->columnScale();
     4763          const int * COIN_RESTRICT column = rowCopy->getIndices();
     4764          const CoinBigIndex * COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
     4765          double * COIN_RESTRICT element = rowCopy->getMutableElements();
     4766          const double * COIN_RESTRICT rowScale = model->rowScale();
     4767          const double * COIN_RESTRICT columnScale = model->columnScale();
    47684768          // scale row copy
    47694769          for (int iRow = 0; iRow < numberRows; iRow++) {
    47704770               CoinBigIndex j;
    47714771               double scale = rowScale[iRow];
    4772                double * elementsInThisRow = element + rowStart[iRow];
    4773                const int * columnsInThisRow = column + rowStart[iRow];
     4772               double * COIN_RESTRICT elementsInThisRow = element + rowStart[iRow];
     4773               const int * COIN_RESTRICT columnsInThisRow = column + rowStart[iRow];
    47744774               int number = rowStart[iRow+1] - rowStart[iRow];
    47754775               assert (number <= numberColumns);
     
    47934793     int numberColumns = matrix_->getNumCols();
    47944794     ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
    4795      const int * row = copy->getIndices();
    4796      const CoinBigIndex * columnStart = copy->getVectorStarts();
    4797      const int * length = copy->getVectorLengths();
    4798      double * element = copy->getMutableElements();
    4799      const double * rowScale = model->rowScale();
    4800      const double * columnScale = model->columnScale();
     4795     const int * COIN_RESTRICT row = copy->getIndices();
     4796     const CoinBigIndex * COIN_RESTRICT columnStart = copy->getVectorStarts();
     4797     const int * COIN_RESTRICT length = copy->getVectorLengths();
     4798     double * COIN_RESTRICT element = copy->getMutableElements();
     4799     const double * COIN_RESTRICT rowScale = model->rowScale();
     4800     const double * COIN_RESTRICT columnScale = model->columnScale();
    48014801     // scale column copy
    48024802     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    48034803          CoinBigIndex j;
    48044804          double scale = columnScale[iColumn];
    4805           double * elementsInThisColumn = element + columnStart[iColumn];
    4806           const int * rowsInThisColumn = row + columnStart[iColumn];
     4805          double * COIN_RESTRICT elementsInThisColumn = element + columnStart[iColumn];
     4806          const int * COIN_RESTRICT rowsInThisColumn = row + columnStart[iColumn];
    48074807          int number = length[iColumn];
    48084808          assert (number <= numberRows);
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r1836 r2078  
    1010
    1111#include "ClpMatrixBase.hpp"
    12 
    13 // Compilers can produce better code if they know about __restrict
    14 #ifndef COIN_RESTRICT
    15 #ifdef COIN_USE_RESTRICT
    16 #define COIN_RESTRICT __restrict
    17 #else
    18 #define COIN_RESTRICT
    19 #endif
    20 #endif
    2112
    2213/** This implements CoinPackedMatrix as derived from ClpMatrixBase.
  • trunk/Clp/src/ClpPlusMinusOneMatrix.cpp

    r2030 r2078  
    1717#include "ClpPlusMinusOneMatrix.hpp"
    1818#include "ClpMessage.hpp"
    19 
     19#ifdef CLP_PLUS_ONE_MATRIX
     20static int oneitcount[13]={0,0,0,0,0,0,0,0,0,0,0,0,0};
     21static void oneit(int i)
     22{
     23  if (!oneitcount[i])
     24    printf("Plus ones for call %d\n",i);
     25  oneitcount[i]++;
     26  oneitcount[12]++;
     27  if ((oneitcount[12]%1000)==0) {
     28    printf("Plus counts");
     29    for (int j=0;j<12;j++)
     30      printf(" %d",oneitcount[j]);
     31    printf("\n");
     32  }   
     33}
     34#endif
    2035//#############################################################################
    2136// Constructors / Destructor / Assignment
     
    3651     numberRows_ = 0;
    3752     numberColumns_ = 0;
     53#ifdef CLP_PLUS_ONE_MATRIX
     54     otherFlags_ = 0;
     55#endif
    3856     columnOrdered_ = true;
    3957}
     
    5270     numberRows_ = rhs.numberRows_;
    5371     numberColumns_ = rhs.numberColumns_;
     72#ifdef CLP_PLUS_ONE_MATRIX
     73     otherFlags_ = rhs.otherFlags_;
     74#endif
    5475     columnOrdered_ = rhs.columnOrdered_;
    5576     if (numberColumns_) {
     
    82103     numberColumns_ = numberColumns;
    83104     columnOrdered_ = columnOrdered;
     105#ifdef CLP_PLUS_ONE_MATRIX
     106     otherFlags_ = 0;
     107#endif
    84108     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
    85109     int numberElements = startPositive[numberMajor];
     
    203227          numberColumns_ = rhs.numberColumns_;
    204228          columnOrdered_ = rhs.columnOrdered_;
     229#ifdef CLP_PLUS_ONE_MATRIX
     230          otherFlags_ = rhs.otherFlags_;
     231#endif
    205232          if (numberColumns_) {
    206233               CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
     
    248275     numberColumns_ = 0;
    249276     columnOrdered_ = rhs.columnOrdered_;
     277#ifdef CLP_PLUS_ONE_MATRIX
     278     otherFlags_ = rhs.otherFlags_;
     279#endif
    250280     if (numberRows <= 0 || numberColumns <= 0) {
    251281          startPositive_ = new CoinBigIndex [1];
     
    427457     return newCopy;
    428458}
     459static bool doPlusOnes=true;
    429460//unscaled versions
    430461void
    431462ClpPlusMinusOneMatrix::times(double scalar,
    432                              const double * x, double * y) const
     463                             const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
    433464{
    434465     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
     
    436467     CoinBigIndex j;
    437468     assert (columnOrdered_);
    438      for (i = 0; i < numberMajor; i++) {
     469#ifdef CLP_PLUS_ONE_MATRIX
     470     if ((otherFlags_&1)==0||!doPlusOnes) {
     471#endif
     472       for (i = 0; i < numberMajor; i++) {
    439473          double value = scalar * x[i];
    440474          if (value) {
     
    448482               }
    449483          }
    450      }
     484       }
     485#ifdef CLP_PLUS_ONE_MATRIX
     486     } else {
     487       // plus one
     488       oneit(0);
     489       for (i = 0; i < numberMajor; i++) {
     490          double value = scalar * x[i];
     491          if (value) {
     492               for (j = startPositive_[i]; j < startPositive_[i+1]; j++) {
     493                    int iRow = indices_[j];
     494                    y[iRow] += value;
     495               }
     496          }
     497       }
     498     }
     499#endif
    451500}
    452501void
    453502ClpPlusMinusOneMatrix::transposeTimes(double scalar,
    454                                       const double * x, double * y) const
     503                                      const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
    455504{
    456505     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
     
    458507     CoinBigIndex j = 0;
    459508     assert (columnOrdered_);
    460      for (i = 0; i < numberMajor; i++) {
     509#ifdef CLP_PLUS_ONE_MATRIX
     510     if ((otherFlags_&1)==0||!doPlusOnes) {
     511#endif
     512       for (i = 0; i < numberMajor; i++) {
    461513          double value = 0.0;
    462514          for (; j < startNegative_[i]; j++) {
     
    469521          }
    470522          y[i] += scalar * value;
    471      }
     523       }
     524#ifdef CLP_PLUS_ONE_MATRIX
     525     } else {
     526       // plus one
     527       oneit(1);
     528       for (i = 0; i < numberMajor; i++) {
     529          double value = 0.0;
     530          for (; j < startPositive_[i+1]; j++) {
     531               int iRow = indices_[j];
     532               value += x[iRow];
     533          }
     534          y[i] += scalar * value;
     535       }
     536     }
     537#endif
    472538}
    473539void
    474540ClpPlusMinusOneMatrix::times(double scalar,
    475                              const double * x, double * y,
     541                             const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
    476542                             const double * /*rowScale*/,
    477543                             const double * /*columnScale*/) const
     
    482548void
    483549ClpPlusMinusOneMatrix::transposeTimes( double scalar,
    484                                        const double * x, double * y,
     550                                       const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
    485551                                       const double * /*rowScale*/,
    486552                                       const double * /*columnScale*/,
     
    500566     // we know it is not scaled
    501567     columnArray->clear();
    502      double * pi = rowArray->denseVector();
     568     double * COIN_RESTRICT pi = rowArray->denseVector();
    503569     int numberNonZero = 0;
    504      int * index = columnArray->getIndices();
    505      double * array = columnArray->denseVector();
     570     int * COIN_RESTRICT index = columnArray->getIndices();
     571     double * COIN_RESTRICT array = columnArray->denseVector();
    506572     int numberInRowArray = rowArray->getNumElements();
    507573     // maybe I need one in OsiSimplex
     
    539605               // need to expand pi into y
    540606               assert(y->capacity() >= numberRows);
    541                double * piOld = pi;
     607               double * COIN_RESTRICT piOld = pi;
    542608               pi = y->denseVector();
    543                const int * whichRow = rowArray->getIndices();
     609               const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    544610               int i;
    545611               // modify pi so can collapse to one loop
     
    548614                    pi[iRow] = scalar * piOld[i];
    549615               }
    550                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     616#ifdef CLP_PLUS_ONE_MATRIX
     617               if ((otherFlags_&1)==0||!doPlusOnes) {
     618#endif
     619                 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    551620                    double value = 0.0;
    552621                    for (; j < startNegative_[iColumn]; j++) {
     
    562631                         index[numberNonZero++] = iColumn;
    563632                    }
    564                }
     633                 }
     634#ifdef CLP_PLUS_ONE_MATRIX
     635               } else {
     636                 // plus one
     637                 oneit(2);
     638                 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     639                    double value = 0.0;
     640                    for (; j < startPositive_[iColumn+1]; j++) {
     641                         int iRow = indices_[j];
     642                         value += pi[iRow];
     643                    }
     644                    if (fabs(value) > zeroTolerance) {
     645                         array[numberNonZero] = value;
     646                         index[numberNonZero++] = iColumn;
     647                    }
     648                 }
     649               }
     650#endif
    565651               for (i = 0; i < numberInRowArray; i++) {
    566652                    int iRow = whichRow[i];
     
    568654               }
    569655          } else {
    570                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     656                 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    571657                    double value = 0.0;
    572658                    for (; j < startNegative_[iColumn]; j++) {
     
    583669                         array[iColumn] = value;
    584670                    }
    585               }
     671                }
    586672          }
    587673          columnArray->setNumElements(numberNonZero);
     
    600686{
    601687     columnArray->clear();
    602      double * pi = rowArray->denseVector();
     688     double * COIN_RESTRICT pi = rowArray->denseVector();
    603689     int numberNonZero = 0;
    604      int * index = columnArray->getIndices();
    605      double * array = columnArray->denseVector();
     690     int * COIN_RESTRICT index = columnArray->getIndices();
     691     double * COIN_RESTRICT array = columnArray->denseVector();
    606692     int numberInRowArray = rowArray->getNumElements();
    607693     // maybe I need one in OsiSimplex
    608694     double zeroTolerance = model->zeroTolerance();
    609      const int * column = indices_;
    610      const CoinBigIndex * startPositive = startPositive_;
    611      const CoinBigIndex * startNegative = startNegative_;
    612      const int * whichRow = rowArray->getIndices();
     695     const int * COIN_RESTRICT column = indices_;
     696     const CoinBigIndex * COIN_RESTRICT startPositive = startPositive_;
     697     const CoinBigIndex * COIN_RESTRICT startNegative = startNegative_;
     698     const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    613699     bool packed = rowArray->packedMode();
    614700     if (numberInRowArray > 2) {
    615701          // do by rows
    616702          int iRow;
    617           double * markVector = y->denseVector(); // probably empty .. but
     703          double * COIN_RESTRICT markVector = y->denseVector(); // probably empty .. but
    618704          int numberOriginal = 0;
    619705          int i;
    620706          if (packed) {
     707               int numberCovered=0;
     708               int numberColumns = getNumCols();
     709               bool sparse=true;
     710               int target=1*numberColumns;
     711               for (int i = 0; i < numberInRowArray; i++) {
     712                 int iRow = whichRow[i];
     713                 numberCovered += startPositive[iRow+1] - startPositive[iRow];
     714                 if (numberCovered>target) {
     715                   sparse=false;
     716                   break;
     717                 }
     718               }
    621719               numberNonZero = 0;
    622                // and set up mark as char array
    623                char * marked = reinterpret_cast<char *> (index + columnArray->capacity());
    624                double * array2 = y->denseVector();
     720               if (sparse) {
     721                 // and set up mark as char array
     722                 char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + columnArray->capacity());
     723                 double * COIN_RESTRICT array2 = y->denseVector();
    625724#ifdef CLP_DEBUG
    626                int numberColumns = model->numberColumns();
    627                for (i = 0; i < numberColumns; i++) {
     725                 int numberColumns = model->numberColumns();
     726                 for (int i = 0; i < numberColumns; i++) {
    628727                    assert(!marked[i]);
    629728                    assert(!array2[i]);
    630                }
    631 #endif
    632                for (i = 0; i < numberInRowArray; i++) {
    633                     iRow = whichRow[i];
    634                     double value = pi[i] * scalar;
    635                     CoinBigIndex j;
    636                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
     729                 }
     730#endif
     731#ifdef CLP_PLUS_ONE_MATRIX
     732                 if ((otherFlags_&1)==0||!doPlusOnes) {
     733#endif
     734                   for (int i = 0; i < numberInRowArray; i++) {
     735                     iRow = whichRow[i];
     736                     double value = pi[i] * scalar;
     737                     CoinBigIndex j;
     738                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
    637739                         int iColumn = column[j];
    638740                         if (!marked[iColumn]) {
     
    641743                         }
    642744                         array2[iColumn] += value;
    643                     }
    644                     for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
     745                      }
     746                      for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
    645747                         int iColumn = column[j];
    646748                         if (!marked[iColumn]) {
     
    649751                         }
    650752                         array2[iColumn] -= value;
    651                     }
    652                }
    653                // get rid of tiny values and zero out marked
    654                numberOriginal = numberNonZero;
    655                numberNonZero = 0;
    656                for (i = 0; i < numberOriginal; i++) {
     753                      }
     754                    }
     755#ifdef CLP_PLUS_ONE_MATRIX
     756                 } else {
     757                    // plus one
     758                    oneit(4);
     759                   for (int i = 0; i < numberInRowArray; i++) {
     760                     iRow = whichRow[i];
     761                     double value = pi[i] * scalar;
     762                     CoinBigIndex j;
     763                     for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) {
     764                         int iColumn = column[j];
     765                         if (!marked[iColumn]) {
     766                              marked[iColumn] = 1;
     767                              index[numberNonZero++] = iColumn;
     768                         }
     769                         array2[iColumn] += value;
     770                      }
     771                    }
     772                }
     773#endif
     774                // get rid of tiny values and zero out marked
     775                 numberOriginal = numberNonZero;
     776                 numberNonZero = 0;
     777                 for (int i = 0; i < numberOriginal; i++) {
    657778                    int iColumn = index[i];
    658779                    if (marked[iColumn]) {
     
    665786                         }
    666787                    }
    667                }
     788                 }
     789               } else {
     790                 // not sparse
     791#ifdef CLP_PLUS_ONE_MATRIX
     792                 if ((otherFlags_&1)==0||!doPlusOnes) {
     793#endif
     794                   for (int i = 0; i < numberInRowArray; i++) {
     795                     iRow = whichRow[i];
     796                     double value = pi[i] * scalar;
     797                     CoinBigIndex j;
     798                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
     799                         int iColumn = column[j];
     800                         array[iColumn] += value;
     801                     }
     802                     for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
     803                         int iColumn = column[j];
     804                         array[iColumn] -= value;
     805                     }
     806                   }
     807#ifdef CLP_PLUS_ONE_MATRIX
     808                } else {
     809                  // plus one
     810                  oneit(5);
     811                   for (int i = 0; i < numberInRowArray; i++) {
     812                     iRow = whichRow[i];
     813                     double value = pi[i] * scalar;
     814                     CoinBigIndex j;
     815                     for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) {
     816                         int iColumn = column[j];
     817                         array[iColumn] += value;
     818                     }
     819                   }
     820                 }
     821#endif
     822                 // get rid of tiny values and count
     823                 for (int i = 0; i < numberColumns; i++) {
     824                   double value = array[i];
     825                   if (value) {
     826                     array[i] = 0.0;
     827                     if (fabs(value) > zeroTolerance) {
     828                       array[numberNonZero] = value;
     829                       index[numberNonZero++] = i;
     830                     }
     831                   }
     832                 }
     833               }
    668834          } else {
    669835               numberNonZero = 0;
    670836               // and set up mark as char array
    671                char * marked = reinterpret_cast<char *> (markVector);
     837               char * COIN_RESTRICT marked = reinterpret_cast<char *> (markVector);
    672838               for (i = 0; i < numberOriginal; i++) {
    673839                    int iColumn = index[i];
     
    675841               }
    676842               for (i = 0; i < numberInRowArray; i++) {
    677                     iRow = whichRow[i];
    678                     double value = pi[iRow] * scalar;
    679                     CoinBigIndex j;
    680                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
    681                          int iColumn = column[j];
    682                          if (!marked[iColumn]) {
    683                               marked[iColumn] = 1;
    684                               index[numberNonZero++] = iColumn;
    685                          }
    686                          array[iColumn] += value;
    687                     }
    688                     for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
    689                          int iColumn = column[j];
    690                          if (!marked[iColumn]) {
    691                               marked[iColumn] = 1;
    692                               index[numberNonZero++] = iColumn;
    693                          }
    694                          array[iColumn] -= value;
    695                     }
    696                }
     843                iRow = whichRow[i];
     844                double value = pi[iRow] * scalar;
     845                CoinBigIndex j;
     846                for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
     847                   int iColumn = column[j];
     848                   if (!marked[iColumn]) {
     849                     marked[iColumn] = 1;
     850                     index[numberNonZero++] = iColumn;
     851                   }
     852                   array[iColumn] += value;
     853                }
     854                for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
     855                   int iColumn = column[j];
     856                   if (!marked[iColumn]) {
     857                     marked[iColumn] = 1;
     858                     index[numberNonZero++] = iColumn;
     859                   }
     860                   array[iColumn] -= value;
     861                }
     862               }
    697863               // get rid of tiny values and zero out marked
    698864               numberOriginal = numberNonZero;
     
    726892               }
    727893               // and set up mark as char array
    728                char * marked = reinterpret_cast<char *> (index + columnArray->capacity());
    729                int * lookup = y->getIndices();
     894               char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + columnArray->capacity());
     895               int * COIN_RESTRICT lookup = y->getIndices();
    730896               double value = pi0 * scalar;
    731                for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
     897               int numberOriginal;
     898#ifdef CLP_PLUS_ONE_MATRIX
     899               if ((otherFlags_&1)==0||!doPlusOnes) {
     900#endif
     901                 for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
    732902                    int iColumn = column[j];
    733903                    array[numberNonZero] = value;
     
    735905                    lookup[iColumn] = numberNonZero;
    736906                    index[numberNonZero++] = iColumn;
    737               }
    738               for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
     907                }
     908                for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
    739909                    int iColumn = column[j];
    740910                    array[numberNonZero] = -value;
     
    742912                    lookup[iColumn] = numberNonZero;
    743913                    index[numberNonZero++] = iColumn;
    744               }
    745                int numberOriginal = numberNonZero;
    746               value = pi1 * scalar;
    747               for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
     914                }
     915                numberOriginal = numberNonZero;
     916                value = pi1 * scalar;
     917                for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
    748918                    int iColumn = column[j];
    749919                    if (marked[iColumn]) {
     
    756926                         }
    757927                    }
    758               }
    759               for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
     928                }
     929                for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
    760930                    int iColumn = column[j];
    761931                    if (marked[iColumn]) {
     
    768938                         }
    769939                    }
    770                }
     940                 }
     941#ifdef CLP_PLUS_ONE_MATRIX
     942               } else {
     943                 // plus one
     944                 oneit(7);
     945                 for (j = startPositive[iRow0]; j < startPositive[iRow0+1]; j++) {
     946                    int iColumn = column[j];
     947                    array[numberNonZero] = value;
     948                    marked[iColumn] = 1;
     949                    lookup[iColumn] = numberNonZero;
     950                    index[numberNonZero++] = iColumn;
     951                 }
     952                 numberOriginal = numberNonZero;
     953                 value = pi1 * scalar;
     954                 for (j = startPositive[iRow1]; j < startPositive[iRow1+1]; j++) {
     955                    int iColumn = column[j];
     956                    if (marked[iColumn]) {
     957                         int iLookup = lookup[iColumn];
     958                         array[iLookup] += value;
     959                    } else {
     960                         if (fabs(value) > zeroTolerance) {
     961                              array[numberNonZero] = value;
     962                              index[numberNonZero++] = iColumn;
     963                         }
     964                    }
     965                 }
     966               }
     967#endif
    771968               // get rid of tiny values and zero out marked
    772969               int nDelete = 0;
     
    8031000               value = pi[iRow0] * scalar;
    8041001               CoinBigIndex j;
    805                for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
    806                     int iColumn = column[j];
    807                     index[numberNonZero++] = iColumn;
    808                     array[iColumn] = value;
    809                }
    810                for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
    811                     int iColumn = column[j];
    812                     index[numberNonZero++] = iColumn;
    813                     array[iColumn] = -value;
    814                }
    815                value = pi[iRow1] * scalar;
    816                for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
    817                     int iColumn = column[j];
    818                     double value2 = array[iColumn];
    819                     if (value2) {
    820                          value2 += value;
    821                     } else {
    822                          value2 = value;
    823                          index[numberNonZero++] = iColumn;
    824                     }
    825                     array[iColumn] = value2;
    826                }
    827                for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
    828                     int iColumn = column[j];
    829                     double value2 = array[iColumn];
    830                     if (value2) {
    831                          value2 -= value;
    832                     } else {
    833                          value2 = -value;
    834                          index[numberNonZero++] = iColumn;
    835                     }
    836                     array[iColumn] = value2;
    837                }
     1002               for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
     1003                int iColumn = column[j];
     1004                index[numberNonZero++] = iColumn;
     1005                array[iColumn] = value;
     1006               }
     1007               for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
     1008                int iColumn = column[j];
     1009                index[numberNonZero++] = iColumn;
     1010                array[iColumn] = -value;
     1011               }
     1012               value = pi[iRow1] * scalar;
     1013               for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
     1014                int iColumn = column[j];
     1015                double value2 = array[iColumn];
     1016                if (value2) {
     1017                   value2 += value;
     1018                } else {
     1019                   value2 = value;
     1020                   index[numberNonZero++] = iColumn;
     1021                }
     1022                array[iColumn] = value2;
     1023               }
     1024               for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
     1025                int iColumn = column[j];
     1026                double value2 = array[iColumn];
     1027                if (value2) {
     1028                   value2 -= value;
     1029                } else {
     1030                   value2 = -value;
     1031                   index[numberNonZero++] = iColumn;
     1032                }
     1033                array[iColumn] = value2;
     1034               }
    8381035               // get rid of tiny values and zero out marked
    8391036               numberOriginal = numberNonZero;
     
    8581055               value = pi[0] * scalar;
    8591056               if (fabs(value) > zeroTolerance) {
     1057#ifdef CLP_PLUS_ONE_MATRIX
     1058                 if ((otherFlags_&1)==0||!doPlusOnes) {
     1059#endif
    8601060                    for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
    8611061                         int iColumn = column[j];
     
    8681068                         index[numberNonZero++] = iColumn;
    8691069                    }
     1070#ifdef CLP_PLUS_ONE_MATRIX
     1071                 } else {
     1072                   // plus one
     1073                   oneit(9);
     1074                    for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) {
     1075                         int iColumn = column[j];
     1076                         array[numberNonZero] = value;
     1077                         index[numberNonZero++] = iColumn;
     1078                    }
     1079                 }
     1080#endif
    8701081               }
    8711082          } else {
    8721083               value = pi[iRow] * scalar;
    8731084               if (fabs(value) > zeroTolerance) {
    874                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
     1085                  for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
    8751086                         int iColumn = column[j];
    8761087                         array[iColumn] = value;
     
    8991110{
    9001111     columnArray->clear();
    901      double * pi = rowArray->denseVector();
    902      double * array = columnArray->denseVector();
     1112     double * COIN_RESTRICT pi = rowArray->denseVector();
     1113     double * COIN_RESTRICT array = columnArray->denseVector();
    9031114     int jColumn;
    9041115     int numberToDo = y->getNumElements();
    905      const int * which = y->getIndices();
     1116     const int * COIN_RESTRICT which = y->getIndices();
    9061117     assert (!rowArray->packedMode());
    9071118     columnArray->setPacked();
    908      for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     1119#ifdef CLP_PLUS_ONE_MATRIX
     1120     if ((otherFlags_&1)==0||!doPlusOnes) {
     1121#endif
     1122       for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    9091123          int iColumn = which[jColumn];
    9101124          double value = 0.0;
     
    9191133          }
    9201134          array[jColumn] = value;
    921      }
     1135       }
     1136#ifdef CLP_PLUS_ONE_MATRIX
     1137     } else {
     1138       // plus one
     1139       oneit(11);
     1140       for (jColumn = 0; jColumn < numberToDo; jColumn++) {
     1141          int iColumn = which[jColumn];
     1142          double value = 0.0;
     1143          CoinBigIndex j = startPositive_[iColumn];
     1144          for (; j < startPositive_[iColumn+1]; j++) {
     1145               int iRow = indices_[j];
     1146               value += pi[iRow];
     1147          }
     1148          array[jColumn] = value;
     1149       }
     1150     }
     1151#endif
    9221152}
    9231153/// returns number of elements in column part of basis,
     
    9901220                                    int iColumn) const
    9911221{
    992      int * index = rowArray->getIndices();
    993      double * array = rowArray->denseVector();
     1222     int * COIN_RESTRICT index = rowArray->getIndices();
     1223     double * COIN_RESTRICT array = rowArray->denseVector();
    9941224     int number = 0;
    9951225     CoinBigIndex j = startPositive_[iColumn];
     
    12491479{
    12501480     columnOrdered_ = columnOrdered;
     1481#ifdef CLP_PLUS_ONE_MATRIX
     1482     otherFlags_ = 0;
     1483#endif
    12511484     startPositive_ = startPositive;
    12521485     startNegative_ = startNegative;
     
    12671500     CoinBigIndex last = -1;
    12681501     int bad = 0;
     1502     // say if all +1
     1503#ifdef CLP_PLUS_ONE_MATRIX
     1504     otherFlags_ = 1;
     1505#endif
    12691506     for (int i = 0; i < number; i++) {
    12701507          if(startPositive_[i] < last)
     
    12721509          else
    12731510               last = startPositive_[i];
     1511#ifdef CLP_PLUS_ONE_MATRIX
     1512          if (startNegative_[i] < startPositive_[i+1])
     1513            otherFlags_ &= ~1; // not all +1
     1514#endif
    12741515          if(startNegative_[i] < last)
    12751516               bad++;
     
    13241565     CoinBigIndex size = 0;
    13251566     int numberBad = 0;
     1567     // say if all +1
     1568#ifdef CLP_PLUS_ONE_MATRIX
     1569     otherFlags_ = 0;
     1570#endif
    13261571     for (iColumn = 0; iColumn < number; iColumn++) {
    13271572          int n = columns[iColumn]->getNumElements();
     
    13881633     CoinBigIndex size = 0;
    13891634     int numberBad = 0;
     1635     // say if all +1
     1636#ifdef CLP_PLUS_ONE_MATRIX
     1637     otherFlags_ = 0;
     1638#endif
    13901639     for (iRow = 0; iRow < number; iRow++) {
    13911640          int n = rows[iRow]->getNumElements();
     
    15091758     CoinBigIndex j;
    15101759     double tolerance = model->currentDualTolerance();
    1511      double * reducedCost = model->djRegion();
    1512      const double * duals = model->dualRowSolution();
    1513      const double * cost = model->costRegion();
     1760     double * COIN_RESTRICT reducedCost = model->djRegion();
     1761     const double * COIN_RESTRICT duals = model->dualRowSolution();
     1762     const double * COIN_RESTRICT cost = model->costRegion();
    15141763     double bestDj;
    15151764     if (bestSequence >= 0)
     
    16411890     matrix_ = NULL;
    16421891     lengths_ = NULL;
     1892     // say if all +1
     1893#ifdef CLP_PLUS_ONE_MATRIX
     1894     otherFlags_ = 0;
     1895#endif
    16431896}
    16441897/* Returns true if can combine transposeTimes and subsetTransposeTimes
     
    16811934                                       double referenceIn, double devex,
    16821935                                       // Array for exact devex to say what is in reference framework
    1683                                        unsigned int * reference,
    1684                                        double * weights, double scaleFactor)
     1936                                       unsigned int * COIN_RESTRICT reference,
     1937                                       double * COIN_RESTRICT weights, double scaleFactor)
    16851938{
    16861939     // put row of tableau in dj1
    1687      double * pi = pi1->denseVector();
     1940     double * COIN_RESTRICT pi = pi1->denseVector();
    16881941     int numberNonZero = 0;
    1689      int * index = dj1->getIndices();
    1690      double * array = dj1->denseVector();
     1942     int * COIN_RESTRICT index = dj1->getIndices();
     1943     double * COIN_RESTRICT array = dj1->denseVector();
    16911944     int numberInRowArray = pi1->getNumElements();
    16921945     double zeroTolerance = model->zeroTolerance();
     
    16951948     int iColumn;
    16961949     assert (!spare->getNumElements());
    1697      double * piWeight = pi2->denseVector();
     1950     double * COIN_RESTRICT piWeight = pi2->denseVector();
    16981951     assert (!pi2->packedMode());
    16991952     bool killDjs = (scaleFactor == 0.0);
     
    17041957          // need to expand pi into y
    17051958          assert(spare->capacity() >= model->numberRows());
    1706           double * piOld = pi;
     1959          double * COIN_RESTRICT piOld = pi;
    17071960          pi = spare->denseVector();
    1708           const int * whichRow = pi1->getIndices();
     1961          const int * COIN_RESTRICT whichRow = pi1->getIndices();
    17091962          int i;
    17101963          // modify pi so can collapse to one loop
     
    18262079                                    double referenceIn, double devex,
    18272080                                    // Array for exact devex to say what is in reference framework
    1828                                     unsigned int * reference,
    1829                                     double * weights, double scaleFactor)
     2081                                    unsigned int * COIN_RESTRICT reference,
     2082                                    double * COIN_RESTRICT weights, double scaleFactor)
    18302083{
    18312084     int number = dj1->getNumElements();
    1832      const int * index = dj1->getIndices();
    1833      double * array = dj1->denseVector();
     2085     const int * COIN_RESTRICT index = dj1->getIndices();
     2086     double * COIN_RESTRICT array = dj1->denseVector();
    18342087     assert( dj1->packedMode());
    18352088
    1836      double * piWeight = pi2->denseVector();
     2089     double * COIN_RESTRICT piWeight = pi2->denseVector();
    18372090     bool killDjs = (scaleFactor == 0.0);
    18382091     if (!scaleFactor)
     
    19332186{
    19342187     int numberErrors = 0;
     2188     // say if all +1
     2189#ifdef CLP_PLUS_ONE_MATRIX
     2190     otherFlags_ = 0;
     2191#endif
    19352192     // make into CoinPackedVector
    19362193     CoinPackedVectorBase ** vectors =
  • trunk/Clp/src/ClpPlusMinusOneMatrix.hpp

    r1665 r2078  
    267267     mutable int * lengths_;
    268268     /// Start of +1's for each
    269      CoinBigIndex * startPositive_;
     269     CoinBigIndex * COIN_RESTRICT startPositive_;
    270270     /// Start of -1's for each
    271      CoinBigIndex * startNegative_;
     271     CoinBigIndex * COIN_RESTRICT startNegative_;
    272272     /// Data -1, then +1 rows in pairs (row==-1 if one entry)
    273      int * indices_;
     273     int * COIN_RESTRICT indices_;
    274274     /// Number of rows
    275275     int numberRows_;
    276276     /// Number of columns
    277277     int numberColumns_;
     278#ifdef CLP_PLUS_ONE_MATRIX
     279     /** Other flags (could have columnOrdered_?)
     280         1 bit - says just +1
     281     */
     282     mutable int otherFlags_;
     283#endif
    278284     /// True if column ordered
    279285     bool columnOrdered_;
  • trunk/Clp/src/ClpSimplex.cpp

    r2077 r2078  
    534534          largeValue_ = value;
    535535}
     536double minimumPrimalToleranceZZ=0.0;
     537#define CLP_INFEAS_SAVE 5
     538double averageInfeasZZ[CLP_INFEAS_SAVE];
    536539int
    537540ClpSimplex::gutsOfSolution ( double * givenDuals,
     
    759762     }
    760763     int returnCode=0;
     764     bool notChanged=true;
    761765     if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_<0 || factorization_->pivotTolerance()<0.9899999999) &&
    762766         (oldLargestDualError||oldLargestPrimalError)) {
     
    780784         else if (pivotTolerance<0.98999999)
    781785           factorization_->pivotTolerance(CoinMin(0.99,pivotTolerance*factor));
     786         notChanged=pivotTolerance==factorization_->pivotTolerance();
    782787#ifdef CLP_USEFUL_PRINTOUT
    783788         if (pivotTolerance<0.9899999) {
     
    797802       }
    798803     }
     804     if (progress_.iterationNumber_[0]>0&&
     805         progress_.iterationNumber_[CLP_PROGRESS-1]
     806         -progress_.iterationNumber_[0]<CLP_PROGRESS*3&&
     807         factorization_->pivotTolerance()<0.25&&notChanged) {
     808       double pivotTolerance = factorization_->pivotTolerance();
     809       factorization_->pivotTolerance(pivotTolerance*1.5);
     810#ifdef CLP_USEFUL_PRINTOUT
     811       printf("Changing pivot tolerance from %g to %g - inverting too often\n",
     812              pivotTolerance,factorization_->pivotTolerance());
     813#endif
     814     }
    799815     // Switch off false values pass indicator
    800816     if (!valuesPass && algorithm_ > 0)
     
    805821            numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,sumOfRelaxedPrimalInfeasibilities_,
    806822            numberDualInfeasibilities_,sumDualInfeasibilities_,sumOfRelaxedDualInfeasibilities_);
     823     if ((moreSpecialOptions_&8388608)!=0) {
     824       if (algorithm_<0) {
     825         bool doneFiddling=false;
     826         // Optimization may make exact test iffy
     827         double testTolerance=minimumPrimalToleranceZZ+1.0e-15;
     828         while (!doneFiddling) {
     829           doneFiddling=true;
     830           while( !sumOfRelaxedPrimalInfeasibilities_&&
     831                  primalTolerance_>testTolerance) {
     832             // feasible - adjust tolerance
     833             double saveTolerance=primalTolerance_;
     834             primalTolerance_=CoinMax(0.25*primalTolerance_,
     835                                      minimumPrimalToleranceZZ);
     836             printf("Resetting primal tolerance from %g to %g\n",
     837                    saveTolerance,primalTolerance_);
     838             dblParam_[ClpPrimalTolerance]=primalTolerance_;
     839             moreSpecialOptions_ &= ~8388608;
     840             // redo with switch off
     841             returnCode=gutsOfSolution ( givenDuals,givenPrimals,valuesPass);
     842           }
     843           if(primalTolerance_>testTolerance)
     844             moreSpecialOptions_ |= 8388608; // back on
     845           if ((moreSpecialOptions_&8388608)!=0) {
     846             assert( numberPrimalInfeasibilities_);
     847             // average infeasibility
     848             double average=sumPrimalInfeasibilities_/numberPrimalInfeasibilities_;
     849             double minimum=COIN_DBL_MAX;
     850             double averageTotal=average;
     851             bool firstTime=averageInfeasZZ[0]==COIN_DBL_MAX;
     852             for (int i=0;i<CLP_INFEAS_SAVE-1;i++) {
     853               double value = averageInfeasZZ[i+1];
     854               averageTotal+=value;
     855               averageInfeasZZ[i]=value;
     856               minimum=CoinMin(minimum,value);
     857             }
     858             averageInfeasZZ[CLP_INFEAS_SAVE-1]=average;
     859             averageTotal /= CLP_INFEAS_SAVE;
     860             double oldTolerance=primalTolerance_;
     861             if (averageInfeasZZ[0]!=COIN_DBL_MAX) {
     862               if (firstTime) {
     863                 primalTolerance_=CoinMin(0.1,0.1*averageTotal);
     864                 primalTolerance_ = CoinMin(primalTolerance_,average);
     865               } else if (primalTolerance_>0.1*minimum) {
     866                 primalTolerance_=0.1*minimum;
     867               }
     868               primalTolerance_=
     869                 CoinMax(primalTolerance_,minimumPrimalToleranceZZ);
     870             }
     871             if (primalTolerance_!=oldTolerance) {
     872               printf("Changing primal tolerance from %g to %g\n",
     873                      oldTolerance,primalTolerance_);
     874               moreSpecialOptions_ &= ~8388608;
     875               // redo with switch off
     876               returnCode=gutsOfSolution ( givenDuals,givenPrimals,valuesPass);
     877               if(primalTolerance_>testTolerance)
     878                 moreSpecialOptions_ |= 8388608|4194304;
     879               if( !sumOfRelaxedPrimalInfeasibilities_)
     880                 doneFiddling=false; // over done it
     881             }
     882           }
     883         }
     884       }
     885     }
    807886     return returnCode;
    808887}
     
    9361015     }
    9371016#endif
     1017     //printf ("ZZ0 n before %d",number);
    9381018     if (number)
    9391019          factorization_->updateColumn(workSpace, thisVector);
     1020     //printf(" - after %d\n",thisVector->getNumElements());
    9401021     double * work = workSpace->denseVector();
    9411022#ifdef CLP_DEBUG
     
    10151096               thisVector->setNumElements(number);
    10161097               lastError = largestPrimalError_;
     1098               //printf ("ZZ%d n before %d",iRefine+1,number);
    10171099               factorization_->updateColumn(workSpace, thisVector);
     1100               //printf(" - after %d\n",thisVector->getNumElements());
    10181101               multiplier = 1.0 / multiplier;
    10191102               double * previous = lastVector->denseVector();
     
    21832266               }
    21842267          }
     2268#if CLP_FACTORIZATION_NEW_TIMING>1
     2269          factorization_->statsRefactor('M');
     2270#endif
    21852271          return 1;
    21862272     } else if ((factorization_->timeToRefactorize() && !dontInvert)
    21872273                ||invertNow) {
    21882274          //printf("ret after %d pivots\n",factorization_->pivots());
     2275#if CLP_FACTORIZATION_NEW_TIMING>1
     2276          factorization_->statsRefactor('T');
     2277#endif
    21892278          return 1;
    21902279     } else if (forceFactorization_ > 0 &&
     
    21942283          if (forceFactorization_ > factorization_->maximumPivots())
    21952284               forceFactorization_ = -1; //off
     2285#if CLP_FACTORIZATION_NEW_TIMING>1
     2286          factorization_->statsRefactor('F');
     2287#endif
    21962288          return 1;
    21972289     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
     
    44554547       moreSpecialOptions_ &= ~4194304;
    44564548       primalTolerance_=1.0e-7;
    4457        dblParam_[ClpPrimalTolerance]=1.0e-7;
     4549       dblParam_[ClpPrimalTolerance]=primalTolerance_;
    44584550       dualTolerance_=1.0e-7;
    4459        dblParam_[ClpDualTolerance]=1.0e-7;
     4551       dblParam_[ClpDualTolerance]=dualTolerance_;
    44604552     }
    44614553     // ray may be null if in branch and bound
     
    54495541     //#define EXPENSIVE
    54505542#endif
     5543     for (int i=0;i<CLP_INFEAS_SAVE;i++)
     5544       averageInfeasZZ[i]=COIN_DBL_MAX;
    54515545#ifdef EXPENSIVE
    54525546     static int dualCount = 0;
     
    1041910513               frequency = base + cutoff1 / freq0 + (numberRows_ - cutoff1) / freq1;
    1042010514#endif
     10515          //frequency *= 1.05;
    1042110516          setFactorizationFrequency(CoinMin(maximum, frequency));
    1042210517     }
  • trunk/Clp/src/ClpSimplex.hpp

    r2074 r2078  
    12391239         2097152 bit - no primal in fastDual2 if feasible
    12401240         4194304 bit - tolerances have been changed by code
     1241         8388608 bit - tolerances are dynamic (at first)
    12411242     */
    12421243     inline void setMoreSpecialOptions(int value) {
  • trunk/Clp/src/ClpSimplexDual.cpp

    r2074 r2078  
    135135#define NDEBUG_CLP
    136136#endif
     137extern double minimumPrimalToleranceZZ;
     138#define CLP_INFEAS_SAVE 5
     139extern double averageInfeasZZ[CLP_INFEAS_SAVE];
    137140// dual
    138141
     
    594597     if (alphaAccuracy_ != -1.0)
    595598          alphaAccuracy_ = 1.0;
     599     minimumPrimalToleranceZZ=primalTolerance();
    596600     int returnCode = startupSolve(ifValuesPass, saveDuals, startFinishOptions);
    597601     // Save so can see if doing after primal
     
    43454349#endif
    43464350         primalTolerance_=1.0e-6;
     4351         minimumPrimalToleranceZZ=primalTolerance_;
    43474352         dblParam_[ClpPrimalTolerance]=1.0e-6;
    43484353         moreSpecialOptions_ |= 4194304;
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r2077 r2078  
    10301030                             numberThrownOut,__FILE__,__LINE__);
    10311031                      //abort();
     1032                    //if (numberThrownOut)
     1033                    //printf("OUCH! - %d thrownout at %s line %d\n",
     1034                    //       numberThrownOut,__FILE__,__LINE__);
    10321035                    sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
    10331036               }
     
    11821185               progress->initialWeight_ = -1.0;
    11831186          }
    1184           //printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
    1185           //   numberSame,numberDifferent,numberZero,
    1186           //   numberFreeSame,numberFreeDifferent,numberFreeZero);
     1187#ifdef CLP_USEFUL_PRINTOUT
     1188          printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
     1189             numberSame,numberDifferent,numberZero,
     1190             numberFreeSame,numberFreeDifferent,numberFreeZero);
     1191#endif
    11871192          // we want most to be same
    11881193          if (n) {
     1194               std::sort(dj_, dj_ + n);
     1195#ifdef CLP_USEFUL_PRINTOUT
    11891196               double most = 0.95;
    1190                std::sort(dj_, dj_ + n);
    11911197               int which = static_cast<int> ((1.0 - most) * static_cast<double> (n));
    1192                double take = -dj_[which] * infeasibilityCost_;
    1193                //printf("XXXXZ inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take,-dj_[0]*infeasibilityCost_,-dj_[n-1]*infeasibilityCost_);
    1194                take = -dj_[0] * infeasibilityCost_;
    1195                infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);;
    1196                //printf("XXXX increasing weight to %g\n",infeasibilityCost_);
     1198               double take2 = -dj_[which] * infeasibilityCost_;
     1199               printf("XXXX inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take2,-dj_[0]*infeasibilityCost_,-dj_[n-1]*infeasibilityCost_);
     1200#endif
     1201               double take = -dj_[0] * infeasibilityCost_;
     1202               // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);
     1203               infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10);
     1204#ifdef CLP_USEFUL_PRINTOUT
     1205               printf("XXXX changing weight to %g\n",infeasibilityCost_);
     1206#endif
    11971207          }
    11981208          delete [] dj_;
  • trunk/Clp/src/ClpSolve.cpp

    r2074 r2078  
    562562  if (!this->abcState()||!numberRows_||!numberColumns_) {
    563563    //this->readBasis("aaa.bas");
    564     if (!solveType)
     564    if (!solveType) {
    565565      this->dual(0);
    566     else if (solveType==1)
    567       this->primal(startUp ? 1 : 0);
     566    } else if (solveType==1) {
     567      int ifValuesPass=startUp ? 1 : 0;
     568      if (startUp==3)
     569        ifValuesPass=2;
     570      this->primal(ifValuesPass);
     571    }
    568572    //this->writeBasis("a.bas",true);
    569573  } else {
     
    579583    abcModel2->setAbcState(CLP_ABC_WANTED|crashState|(abcModel2->abcState()&15));
    580584    int ifValuesPass=startUp ? 1 : 0;
     585    if (startUp==3)
     586      ifValuesPass=2;
    581587    // temp
    582588    if (fabs(this->primalTolerance()-1.001e-6)<0.999e-9) {
     
    617623        abcModel2->setStateOfProblem(abcModel2->stateOfProblem() | VALUES_PASS);
    618624        Idiot info(*abcModel2);
     625        info.setStrategy(idiotOptions | info.getStrategy());
    619626        info.setStrategy(512 | info.getStrategy());
    620627        // Allow for scaling
     
    697704    } else if (solveType==1) {
    698705      int saveOptions=abcModel2->specialOptions();
    699       if (startUp==2)
    700         abcModel2->setSpecialOptions(8192|saveOptions);
     706      //if (startUp==2)
     707      //abcModel2->setSpecialOptions(8192|saveOptions);
    701708      abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
    702709      abcModel2->setSpecialOptions(saveOptions);
     
    850857     ClpMatrixBase * saveMatrix = NULL;
    851858     ClpObjective * savedObjective = NULL;
     859     int idiotOptions=0;
     860     if(options.getSpecialOption(6))
     861       idiotOptions=options.getExtraInfo(6)*32768;
    852862#ifdef CLP_USEFUL_PRINTOUT
    853863     debugInt[0]=numberRows();
     
    12631273            }
    12641274          }
    1265           if (nNon <= 1 || allOnes) {
     1275          if (nNon <= 1 || allOnes ||
     1276              (options.getExtraInfo(1) > 0 && options.getSpecialOption(1)==2)) {
    12661277#ifdef COIN_DEVELOP
    12671278               printf("Forcing primal\n");
     
    15441555                    // switch off idiot or sprint if any free variable
    15451556                    // switch off sprint if very few with costs
     1557                    // or great variation in cost
    15461558                    int iColumn;
    15471559                    const double * columnLower = model2->columnLower();
     
    15501562                    int nObj=0;
    15511563                    int nFree=0;
     1564                    double smallestObj=COIN_DBL_MAX;
     1565                    double largestObj=0.0;
    15521566                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    15531567                         if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
     
    15551569                         } else if (objective[iColumn]) {
    15561570                           nObj++;
     1571                           smallestObj=CoinMin(smallestObj,objective[iColumn]);
     1572                           largestObj=CoinMax(largestObj,objective[iColumn]);
    15571573                         }
    15581574                    }
    1559                     if (nObj*10<numberColumns)
     1575                    if (nObj*10<numberColumns ||
     1576                        smallestObj*10.0<largestObj)
    15601577                      doSprint=0;
    15611578                    if (nFree)
     
    16331650                         } else {
    16341651                              nPasses = 10 + numberColumns / 1000;
    1635                               nPasses = CoinMin(nPasses, 200);
    16361652                              nPasses = CoinMax(nPasses, 100);
    16371653                              if (!largestGap)
    16381654                                   nPasses *= 2;
     1655                              nPasses = CoinMin(nPasses, 200);
    16391656                         }
    16401657                    }
     
    17601777               nPasses = 0;
    17611778               Idiot info(*model2);
     1779               info.setStrategy(idiotOptions | info.getStrategy());
    17621780               // Get average number of elements per column
    17631781               double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
     
    19381956               int nPasses = 0;
    19391957               Idiot info(*model2);
     1958               info.setStrategy(idiotOptions | info.getStrategy());
    19401959               // Get average number of elements per column
    19411960               double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
     
    20622081                         //info.setStartingWeight(1.0e-1);
    20632082                    }
    2064                     if (tryIt==1)
    2065                       model2->setSpecialOptions(model2->specialOptions()
    2066                                                 |8388608);
     2083                    if (tryIt==1) {
     2084                      idiotOptions |= 262144;
     2085                      info.setStrategy(idiotOptions | info.getStrategy());
     2086                      //model2->setSpecialOptions(model2->specialOptions()
     2087                      //                        |8388608);
     2088                    }
    20672089               }
    20682090               if (doIdiot > 0) {
     
    21712193#ifndef LACI_TRY
    21722194                    //if (doIdiot>0)
    2173 #ifdef ABC_INHERIT
     2195#if 0 //def ABC_INHERIT
    21742196                    if (!model2->abcState())
    21752197#endif
     
    37093731     delete pinfo;
    37103732     moreSpecialOptions_= saveMoreOptions;
     3733#ifdef CLP_USEFUL_PRINTOUT
     3734     debugInt[23]=numberIterations_;
     3735#endif
    37113736     return finalStatus;
    37123737}
  • trunk/Clp/src/ClpSolve.hpp

    r2074 r2078  
    100100         5 - for presolve
    101101                      1 - switch off dual stuff
    102          6 - for detailed printout (initially just presolve)
    103                       1 - presolve statistics
     102         6 - extra switches
     103                     
    104104     */
    105105     void setSpecialOption(int which, int value, int extraInfo = -1);
  • trunk/Clp/src/ClpSolver.cpp

    r2074 r2078  
    213213#define DEFAULT_PRESOLVE_PASSES 20
    214214#else
    215 #define DEFAULT_PRESOLVE_PASSES 5
     215#define DEFAULT_PRESOLVE_PASSES 10
    216216#endif
    217217          int preSolve = DEFAULT_PRESOLVE_PASSES;
     
    942942                                   } else if (type == CLP_PARAM_ACTION_EITHERSIMPLEX) {
    943943                                        method = ClpSolve::automatic;
     944                                        if (doCrash>3) {
     945                                          solveOptions.setSpecialOption(6, 1,doCrash-3);
     946                                          doCrash=0;
     947                                        }
     948                                        if (doIdiot > 0)
     949                                             solveOptions.setSpecialOption(1, 2, doIdiot);
    944950                                   } else {
    945951                                        method = ClpSolve::useBarrier;
     
    970976                                   if (method == ClpSolve::useDual) {
    971977                                        // dual
    972                                         if (doCrash)
     978                                        if (doCrash&&doCrash<4)
    973979                                             solveOptions.setSpecialOption(0, 1, doCrash); // crash
    974980                                        else if (doIdiot)
     
    984990                                             method = ClpSolve::usePrimal;
    985991                                        }
    986                                         if (doCrash) {
     992                                        if (doCrash>3) {
     993                                          solveOptions.setSpecialOption(6, 1,doCrash-3);
     994                                          doCrash=0;
     995                                        }
     996                                        if (doCrash>0) {
    987997                                             solveOptions.setSpecialOption(1, 1, doCrash); // crash
    988998                                        } else if (doSprint > 0) {
     
    29062916#endif
    29072917#ifdef CLP_USEFUL_PRINTOUT
    2908      printf("BENCHMARK %s took %g cpu seconds (%g elapsed) - %d row, %d columns, %d elements - reduced to %d, %d, %d\n",
     2918     printf("BENCHMARK_RESULT %s took %g cpu seconds (%g elapsed - %d iterations) - %d row, %d columns, %d elements - reduced to %d, %d, %d\n",
    29092919              mpsFile.c_str(),CoinCpuTime()-startCpu,CoinGetTimeOfDay()-startElapsed,
     2920            debugInt[23],
    29102921              debugInt[0],debugInt[1],debugInt[2],
    29112922              debugInt[3],debugInt[4],debugInt[5]);
  • trunk/Clp/src/CoinAbcBaseFactorization1.cpp

    r2074 r2078  
    147147#if ABC_SMALL<4
    148148#if ABC_DENSE_CODE>0
    149     denseThreshold_=4;//16; //31;
    150     //denseThreshold_=71;
     149    //denseThreshold_=4;//16; //31;
     150    denseThreshold_=71;
    151151#else
    152152    denseThreshold_=-16;
     
    766766  return status_;
    767767}
     768#define ALWAYS_GIVE_UP 0 //1
    768769#ifdef ABC_USE_FUNCTION_POINTERS
    769770#define DENSE_TRY
     
    776777                         const int * COIN_RESTRICT indexRowU)
    777778{
    778   if (last-first>giveUp) {
     779  if (last-first>giveUp && !ALWAYS_GIVE_UP) {
    779780    int mid=(last+first)>>1;
    780781    cilk_spawn pivotStartup(first,mid,numberInPivotColumn,lengthArea,giveUp,
     
    808809                       CoinFactorizationDouble * COIN_RESTRICT eArea,const CoinFactorizationDouble * COIN_RESTRICT multipliersL)
    809810{
    810   if (last-first>giveUp) {
     811  if (last-first>giveUp && !ALWAYS_GIVE_UP) {
    811812    int mid=(last+first)>>1;
    812813    cilk_spawn pivotWhile(first,mid,numberInPivotColumn,lengthArea,giveUp,
     
    858859                      const CoinFactorizationDouble * COIN_RESTRICT multipliersL, double tolerance)
    859860{
    860   if (last-first>giveUp) {
     861  if (last-first>giveUp && !ALWAYS_GIVE_UP) {
    861862    int mid=(last+first)>>1;
    862863    cilk_spawn pivotSomeAfter(first,mid,numberInPivotColumn,lengthArea,giveUp,
     
    976977                      const CoinFactorizationDouble * COIN_RESTRICT multipliersL, double tolerance)
    977978{
    978   if (last-first>giveUp) {
     979  if (last-first>giveUp && !ALWAYS_GIVE_UP) {
    979980    int mid=(last+first)>>1;
    980981    cilk_spawn pivotSome(first,mid,numberInPivotColumn,lengthArea,giveUp,
     
    24202421      CoinSimplexDouble ratio;
    24212422      if (numberRowsLeft_>2000)
    2422         ratio=4.5;
     2423        ratio=3.5;
    24232424      else if (numberRowsLeft_>800)
    2424         ratio=3.0;//3.5;
     2425        ratio=2.75;//3.5;
    24252426      else if (numberRowsLeft_>256)
    24262427        ratio=2.0;//2.5;
     
    24292430      else
    24302431        ratio=1.5;
     2432      //ratio=200.0;
    24312433      //ratio *=0.75;
    24322434      if ((ratio*leftElements>full&&numberRowsLeft_>abs(denseThreshold_))) {
  • trunk/Clp/src/CoinAbcBaseFactorization4.cpp

    r2024 r2078  
    43664366    if (numberRows_<largeRowsSparse) {
    43674367      sparseThreshold_=CoinMin(numberRows_/mediumRowsDivider,mediumRowsMinCount);
     4368      sparseThreshold_=CoinMin(numberRows_/6,500);
    43684369    } else {
    43694370      sparseThreshold_=CoinMax(largeRowsCount,numberRows_>>3);
     4371      sparseThreshold_=500;
    43704372    }
    43714373#if FACTORIZATION_STATISTICS
  • trunk/Clp/src/IdiSolve.cpp

    r2074 r2078  
    347347     }
    348348     memset(statusSave+numberColumns,0,ncols-numberColumns);
    349      for (int i=0;i<numberColumns;i++) {
    350        if (model_->getColumnStatus(i)==ClpSimplex::isFixed) {
    351          assert (colsol[i]<lower[i]+tolerance||
    352                  colsol[i]>upper[i]-tolerance);
     349     if ((strategy_&131072)==0) {
     350       for (int i=0;i<numberColumns;i++) {
     351         if (model_->getColumnStatus(i)==ClpSimplex::isFixed) {
     352           assert (colsol[i]<lower[i]+tolerance||
     353                   colsol[i]>upper[i]-tolerance);
     354         }
    353355       }
    354356     }
  • trunk/Clp/src/Idiot.cpp

    r2074 r2078  
    7171                      const double * rowLower, const double * rowUpper,
    7272                      const double * cost, const double * element, double fixTolerance,
    73                       double & objValue, double & infValue)
     73                      double & objValue, double & infValue, double & maxInfeasibility)
    7474{
    7575     int n = 0;
     
    111111          objValue = 0.0;
    112112          infValue = 0.0;
     113          maxInfeasibility = 0.0;
    113114          for ( i = 0; i < ncols; i++) {
    114115               if (nextSlack[i] == -1) {
     
    274275                    rowsol[i] = rowValue;
    275276               }
    276                infValue += CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]);
     277               double infeasibility = CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]);
     278               infValue += infeasibility;
     279               maxInfeasibility=CoinMax(maxInfeasibility,infeasibility);
    277280               // just change
    278281               rowsol[i] -= rowSave;
     
    478481#ifndef OSI_IDIOT
    479482     bool fixAfterSome = false; //(model_->specialOptions()&8388608)!=0;
    480      int exitAfter = 1000000; //(model_->specialOptions()&8388608)!=0 ? 50 : 1000000;
     483     int exitAfter = 50; //(model_->specialOptions()&8388608)!=0 ? 50 : 1000000;
    481484     {
    482485        int numberColumns = model_->numberColumns();
     
    717720     }
    718721#endif
     722     if ((strategy_&131072)!=0) {
     723       // going to mess with cost
     724       if (cost==model_->objective())
     725         cost=CoinCopyOfArray(cost,ncols);
     726     }
    719727     for (i = 0; i < ncols; i++) {
    720728          CoinBigIndex j;
     
    757765                           maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
    758766     // update whenUsed_
     767     double maxInfeasibility=COIN_DBL_MAX;
    759768     n = cleanIteration(iteration, ordStart, ordEnd,
    760769                        colsol,  lower,  upper,
    761770                        rowlower, rowupper,
    762                         cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas);
     771                        cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
    763772     if ((strategy_ & 16384) != 0) {
    764773          int * posSlack = whenUsed_ + ncols;
     
    799808               break;
    800809          iteration++;
    801           if (iteration>=exitAfter)
    802             break;
     810          //if (iteration>=exitAfter)
     811          //break;
    803812          checkIteration++;
    804813          if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
     
    826835                             colsol,  lower,  upper,
    827836                             rowlower, rowupper,
    828                              cost, elemXX, fixTolerance, result.objval, result.infeas);
     837                             cost, elemXX, fixTolerance, result.objval, result.infeas, maxInfeasibility);
    829838          if ((strategy_ & 16384) != 0) {
    830839               int * posSlack = whenUsed_ + ncols;
     
    834843               for (i = 0; i < nrows; i++)
    835844                    rowsol[i] += rowsol2[i];
     845          } else {
     846            maxInfeasibility=0.0;
     847            for (i = 0; i < nrows; i++) {
     848              double infeasibility = CoinMax(CoinMax(0.0, rowlower[i] - rowsol[i]), rowsol[i] - rowupper[i]);
     849              maxInfeasibility=CoinMax(maxInfeasibility,infeasibility);
     850            }
    836851          }
    837852          if ((logLevel_ & 1) != 0) {
     
    864879              }
    865880#endif
    866           if (iteration > 50 && n == numberAway ) {
    867             if((result.infeas < 1.0e-4 && majorIterations_<200)||result.infeas<1.0e-8) {
     881          if (iteration > exitAfter ) {
     882            if((result.infeas < 1.0e-4 && majorIterations_<200 && n == numberAway)||result.infeas<1.0e-8||maxInfeasibility<1.0e-8) {
    868883#ifdef CLP_INVESTIGATE
    869884               printf("infeas small %g\n", result.infeas);
    870885#endif
    871                break; // not much happening
     886               if ((strategy_&131072)==0) {
     887                 break; // not much happening
     888               } else {
     889                 int n=0;
     890                 for (int i=0;i<ncols;i++) {
     891                   if (cost[i])
     892                     n++;
     893                 }
     894                 if (n*10<ncols) {
     895                   // fix ones with costs
     896                   exitAfter=100;
     897                   for (int i=0;i<ncols;i++) {
     898                     if (cost[i]||colsol[i]<lower[i]+1.0e-9||
     899                         colsol[i]>upper[i]-1.0e-9) {
     900                       cost[i]=0.0;
     901                       model_->setColumnStatus(i,ClpSimplex::isFixed);
     902                     } else {
     903                       cost[i]=0.0;
     904                       double gap = upper[i]-lower[i];
     905                       if (colsol[i]>upper[i]-0.25*gap) {
     906                         cost[i]=gap/(colsol[i]-upper[i]);
     907                       } else if (colsol[i]<lower[i]+0.25*gap) {
     908                         cost[i]=gap/(colsol[i]-lower[i]);
     909                       }
     910                     }
     911                   }
     912                 }
     913               }
    872914            }
    873915          }
     
    10761118          cost = model_->objective();
    10771119          //rowlower = model_->rowLower();
     1120     } else if (cost != model_->objective()) {
     1121       delete [] cost;
     1122       cost = model_->objective();
    10781123     }
    10791124#endif
     
    11061151     muAtExit_ = mu;
    11071152     // For last iteration make as feasible as possible
    1108      if (oddSlacks && (strategy_&32768)==0)
     1153     if (oddSlacks)
    11091154          strategy_ |= 16384;
    11101155     // not scaled
     
    11121157                        colsol,  lower,  upper,
    11131158                        model_->rowLower(), model_->rowUpper(),
    1114                         cost, element, fixTolerance, lastResult.objval, lastResult.infeas);
     1159                        cost, element, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
    11151160#if 0
    11161161     if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) {
     
    17151760          wantVector = clpMatrixO->wantsSpecialColumnCopy();
    17161761     }
    1717      if ((model_->specialOptions()&8388608)!=0) { // temp test
    1718        justValuesPass=true;
     1762     if ((strategy_&262144)!=0) {
     1763       //justValuesPass=true;
    17191764       assert (addAll<3);
    17201765       assert (presolve);
     1766       //allowInfeasible=true;
     1767     }
     1768     if ((strategy_&32768)!=0)
    17211769       allowInfeasible=true;
    1722      }
    1723      double * saveBounds=NULL;
     1770     if ((strategy_&65536)!=0)
     1771       justValuesPass=true;
     1772     //double * saveBounds=NULL;
    17241773     if (addAll < 3) {
    17251774          ClpPresolve pinfo;
    17261775          if (presolve) {
     1776               double * rhs = new double[nrows];
     1777               double * saveBounds=new double[2*ncols];
     1778               char line[200];
     1779               memcpy(saveBounds,lower,ncols*sizeof(double));
     1780               memcpy(saveBounds+ncols,upper,ncols*sizeof(double));
    17271781               if (allowInfeasible) {
    17281782                    // fix up so will be feasible
    17291783                    const double * dual = model_->dualRowSolution();
    1730                     double * rhs = new double[nrows];
    1731                     double * saveBounds=new double[2*ncols];
    1732                     memcpy(saveBounds,lower,ncols*sizeof(double));
    1733                     memcpy(saveBounds+ncols,upper,ncols*sizeof(double));
    17341784                    for (i = 0; i < nrows; i++)
    17351785                      rhs[i]=fabs(dual[i]);
     
    17371787                    int nSmall=nrows;
    17381788                    int nMedium=nrows;
    1739                     double small=CoinMax(1.0e-4,1.0e-5*rhs[nrows-1]);
    1740                     double medium=CoinMax(1.0e-2,1.0e-3*rhs[nrows-1]);
    1741                     char * marked = new char [nrows];
     1789                    double largest=rhs[nrows-1];
     1790                    double small=CoinMax(1.0e-4,1.0e-5*largest);
     1791                    small = CoinMin(small,1.0e-2);
     1792                    double medium=small*100.0;
    17421793                    double * rowupper = model_->rowUpper();
    17431794                    double * rowlower = model_->rowLower();
     
    17551806                      }
    17561807                    }
    1757                     printf("%d < %g, %d < %g, %d >= %g\n",
    1758                            nSmall,small,nMedium-nSmall,medium,nrows-nMedium,medium);
    1759                     for (i = 0; i < nrows; i++) {
    1760                       if (fabs(dual[i])<medium) {
    1761                         marked[i]=0;
    1762                       } else {
    1763                         marked[i]=1;
     1808                    printf("%d < %g, %d < %g, %d <= %g\n",
     1809                           nSmall,small,nMedium-nSmall,medium,nrows-nMedium,largest);
     1810                    memset(rhs, 0, nrows * sizeof(double));
     1811                    int nFixed=0;
     1812                    for (i = 0; i < ncols; i++) {
     1813                      if (colsol[i]<lower[i]+1.0e-8) {
     1814                        upper[i]=lower[i];
     1815                        colsol[i]=lower[i];
     1816                        nFixed++;
     1817                      } else if (colsol[i]>upper[i]-1.0e-8) {
     1818                        lower[i]=upper[i];
     1819                        colsol[i]=lower[i];
     1820                        nFixed++;
    17641821                      }
    17651822                    }
    1766                     memset(rhs, 0, nrows * sizeof(double));
    17671823                    model_->clpMatrix()->times(1.0, colsol, rhs);
    17681824                    saveRowUpper = CoinCopyOfArray(rowupper, nrows);
    17691825                    saveRowLower = CoinCopyOfArray(rowlower, nrows);
    17701826                    double sum = 0.0;
     1827                    for (i = 0; i < nrows; i++) {
     1828                         if (rhs[i] > rowupper[i]) {
     1829                              sum += rhs[i] - rowupper[i];
     1830                         }
     1831                         if (rhs[i] < rowlower[i]) {
     1832                              sum += rowlower[i] - rhs[i];
     1833                         }
     1834                    }
     1835                    double averageInfeasibility=sum/nrows;
     1836                    double check=CoinMin(1.0e-3,0.1*averageInfeasibility);
    17711837                    int nFixedRows=0;
     1838                    int nFreed=0;
     1839#define MESS_UP 0
    17721840                    for (i = 0; i < nrows; i++) {
    1773                          if (marked[i]&&rowupper[i]>rowlower[i]) {
    1774                            double distance=CoinMin(rowupper[i]-rhs[i],rhs[i]-rowlower[i]);
     1841                         if (rowupper[i]>rowlower[i]+check) {
    17751842                           // look at distance and sign of dual
    1776                            if (rhs[i] > rowupper[i]-1.0e-6) {
    1777                              rhs[i]=rowupper[i];
     1843                           if (dual[i]<-medium&&rowupper[i]-rhs[i]<check) {
     1844                             rowupper[i]=rhs[i];
    17781845                             rowlower[i]=rowupper[i];
    17791846                             nFixedRows++;
    1780                            } else if (rhs[i] < rowlower[i]+1.0e-6) {
    1781                              rhs[i]=rowlower[i];
     1847                           } else if (dual[i]>medium&&rhs[i]-rowlower[i]<check) {
     1848                             rowlower[i]=rhs[i];
    17821849                             rowupper[i]=rowlower[i];
    17831850                             nFixedRows++;
     1851                           } else if (fabs(dual[i])<small&&
     1852                                      rhs[i]-rowlower[i]>check&&
     1853                                      rowupper[i]-rhs[i]>check) {
     1854                             nFreed++;
     1855#if MESS_UP ==1 || MESS_UP ==2
     1856                             rowupper[i]=COIN_DBL_MAX;
     1857                             rowlower[i]=-COIN_DBL_MAX;
     1858#endif
    17841859                           }
    17851860                         }
    17861861                         if (rhs[i] > rowupper[i]) {
    1787                               sum += rhs[i] - rowupper[i];
    17881862                              rowupper[i] = rhs[i];
    17891863                              // maybe make equality
     1864#if MESS_UP ==2 || MESS_UP ==3
     1865                              rowlower[i] = rhs[i];
     1866#endif
    17901867                         }
    17911868                         if (rhs[i] < rowlower[i]) {
    1792                               sum += rowlower[i] - rhs[i];
    17931869                              rowlower[i] = rhs[i];
    17941870                              // maybe make equality
    1795                          }
    1796                     }
     1871#if MESS_UP ==2 || MESS_UP ==3
     1872                              rowupper[i] = rhs[i];
     1873#endif
     1874                         }
     1875                    }
     1876                    sprintf(line,"sum of infeasibilities %g - %d fixed rows, %d fixed columns - might free %d rows",
     1877                           sum,nFixedRows,nFixed,nFreed);
     1878               } else {
     1879                    memset(rhs, 0, nrows * sizeof(double));
    17971880                    int nFixed=0;
    17981881                    for (i = 0; i < ncols; i++) {
    1799                       if (colsol[i]<lower[i]+1.0e-7) {
     1882                      if (colsol[i]<lower[i]+1.0e-8) {
    18001883                        upper[i]=lower[i];
     1884                        colsol[i]=lower[i];
    18011885                        nFixed++;
    1802                       } else if (colsol[i]>upper[i]-1.0e-7) {
     1886                      } else if (colsol[i]>upper[i]-1.0e-8) {
    18031887                        lower[i]=upper[i];
     1888                        colsol[i]=lower[i];
    18041889                        nFixed++;
    18051890                      }
    18061891                    }
    1807                     printf("sum of infeasibilities %g - %d fixed rows, %d fixed columns\n",
    1808                            sum,nFixedRows,nFixed);
    1809                     //COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n", sum));
    1810                     delete [] rhs;
    1811                     delete [] marked;
    1812                }
     1892                    model_->clpMatrix()->times(1.0, colsol, rhs);
     1893                    double sum = 0.0;
     1894                    for (i = 0; i < nrows; i++) {
     1895                         if (rhs[i] > rowupper[i]) {
     1896                              sum += rhs[i] - rowupper[i];
     1897                         }
     1898                         if (rhs[i] < rowlower[i]) {
     1899                              sum += rowlower[i] - rhs[i];
     1900                         }
     1901                    }
     1902                    double averageInfeasibility=sum/nrows;
     1903                    sprintf(line,"sum of infeasibilities %g - average %g, %d fixed columns",
     1904                           sum,averageInfeasibility,nFixed);
     1905               }
     1906               const CoinMessages * messages = model_->messagesPointer();
     1907               model_->messageHandler()->message(CLP_GENERAL, *messages)
     1908                 << line
     1909                 << CoinMessageEol;
     1910               delete [] rhs;
    18131911               saveModel = model_;
    18141912               pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
     
    18191917                 delete [] saveBounds;
    18201918               }
    1821                if (model_&&(saveModel->specialOptions()&8388608)!=0) {
     1919               if (model_&&(strategy_&262144)!=0) {
    18221920                 int nrows = model_->getNumRows();
    18231921                 int ncols = model_->getNumCols();
     
    18511949                 int status=-1;
    18521950                 // If initial is too dense - then all-slack may be better
    1853                  double areaFactor=4.0;
     1951                 double areaFactor=1.0; // was 4.0
    18541952                 const CoinPackedMatrix * matrix = model_->matrix();
    18551953                 while (status) {
     
    18991997          }
    19001998          if (model_) {
     1999               // See if we want to go all way
     2000               int oldSize=2*saveModel->numberRows()+saveModel->numberColumns();
     2001               int newSize=2*model_->numberRows()+model_->numberColumns();
     2002               if (oldSize*2>newSize*3)
     2003                 justValuesPass=false;
    19012004               if (!wantVector) {
    19022005                    //#define TWO_GOES
    19032006#ifdef ABC_INHERIT
    19042007#ifndef TWO_GOES
    1905                     model_->dealWithAbc(1,justValuesPass ? 2 : 1);
     2008                    model_->dealWithAbc(1,justValuesPass ? 3 : 1);
    19062009#else
    19072010                    model_->dealWithAbc(1,1 + 11);
     
    19422045               if (justValuesPass)
    19432046#ifdef ABC_INHERIT
    1944                     model_->dealWithAbc(1,2);
     2047                    model_->dealWithAbc(1,3);
    19452048#else
    19462049                    model_->primal(2);
  • trunk/Clp/src/Idiot.hpp

    r2074 r2078  
    250250                        const double * rowLower, const double * rowUpper,
    251251                        const double * cost, const double * element, double fixTolerance, double & objChange,
    252                         double & infChange);
     252                        double & infChange, double & maxInfeasibility);
    253253private:
    254254     /// Underlying model
     
    287287                  4096 - return best solution (not last found)
    288288                  8192 - always do a presolve in crossover
    289                  16384 - costed slacks found - so whenUsed_ longer */
     289                 16384 - costed slacks found - so whenUsed_ longer
     290                 32768 - experimental 1
     291                 65536 - experimental 2
     292                 131072 - experimental 3
     293                 262144 - just values pass etc */
     294                 
    290295     int lightWeight_; // 0 - normal, 1 lightweight
    291296};
Note: See TracChangeset for help on using the changeset viewer.