Changeset 2146


Ignore:
Timestamp:
Jun 21, 2015 12:23:39 PM (3 years ago)
Author:
forrest
Message:

stuff from stable and start aboca_lite

Location:
trunk/Clp/src
Files:
7 edited

Legend:

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

    r2081 r2146  
    10101010            paction_ = duprow_action::presolve(prob, paction_);
    10111011            printProgress('D',0);
     1012            //paction_ = doubleton_action::presolve(prob, paction_);
     1013            //printProgress('d',0);
     1014            //if (doDependency()) {
     1015              //paction_ = duprow3_action::presolve(prob, paction_);
     1016              //printProgress('Z',0);
     1017            //}
    10121018          }
    10131019          if (doGubrow()) {
     
    10161022               printProgress('E',0);
    10171023          }
     1024          if (ifree) {
     1025            int fill_level=10;
     1026            const CoinPresolveAction * lastAction = NULL;
     1027            int iPass=4;
     1028            while(lastAction!=paction_&&iPass) {
     1029              lastAction=paction_;
     1030              paction_ = implied_free_action::presolve(prob, paction_, fill_level);
     1031              printProgress('l',0);
     1032              iPass--;
     1033            }
     1034          }
    10181035
    10191036          if ((presolveActions_ & 16384) != 0)
     
    13661383#endif
    13671384          }
     1385     }
     1386     if (!prob->status_&&doDependency()) {
     1387       paction_ = duprow3_action::presolve(prob, paction_);
     1388       printProgress('Z',0);
    13681389     }
    13691390     prob->presolveOptions_ &= ~0x10000;
     
    17241745     feasibilityTolerance_(0.0),
    17251746     status_(-1),
     1747     pass_(0),
    17261748     colsToDo_(new int [ncols0_in]),
    17271749     numberColsToDo_(0),
  • trunk/Clp/src/ClpPresolve.hpp

    r2074 r2146  
    152152          else presolveActions_ |= 256;
    153153     }
     154     /// Whether we want to do dependency part of presolve
     155     inline bool doDependency() const {
     156          return (presolveActions_ & 32768) != 0;
     157     }
     158     inline void setDoDependency(bool doDependency) {
     159          if (doDependency) presolveActions_  |= 32768;
     160          else presolveActions_ &= ~32768;
     161     }
    154162     /// Whether we want to do singleton column part of presolve
    155163     inline bool doSingletonColumn() const {
  • trunk/Clp/src/ClpSimplex.cpp

    r2140 r2146  
    1046410464}
    1046510465#endif
     10466/* Clean primal solution
     10467   If you expect solution to only have exact multiples of "exactMultiple" then
     10468   this tries moving solution values to nearest multiple.  If still feasible
     10469   then the solution is replaced.
     10470   
     10471   This is designed for the case where values should be integral, but Clp may
     10472   have values at e.g. 1.0e-13
     10473   Returns 0 if successful, n if n rhs violated
     10474   The dual version will be written if this gets used.
     10475*/
     10476int
     10477ClpSimplex::cleanPrimalSolution(double exactMultiple)
     10478{
     10479  double * tempColumn = new double [numberRows_+numberColumns_];
     10480  double * tempRow = tempColumn+numberColumns_;
     10481  // allow tiny amount of error if not 1.0
     10482  double allowedError=0.0;
     10483  // set up cleaned solution
     10484  if (exactMultiple==1.0) {
     10485    for (int i=0;i<numberColumns_;i++) {
     10486      tempColumn[i] = floor(columnActivity_[i]+0.5);
     10487    }
     10488  } else {
     10489    double reciprocal=1.0/exactMultiple;
     10490    allowedError=1.0e-1*primalTolerance_;
     10491    for (int i=0;i<numberColumns_;i++) {
     10492      double value = floor(columnActivity_[i]*reciprocal+0.5);
     10493      tempColumn[i] = exactMultiple*value;
     10494    }
     10495  }
     10496  // Check bounds
     10497  int nBad=0;
     10498  for (int i=0;i<numberColumns_;i++) {
     10499    double value = tempColumn[i];
     10500    if (value<columnLower_[i]-allowedError||value>columnUpper_[i]+allowedError)
     10501      nBad++;
     10502  }
     10503  memset(tempRow,0,numberRows_*sizeof(double));
     10504  times(-1.0,tempColumn,tempRow);
     10505  for (int i=0;i<numberRows_;i++) {
     10506    double value = tempRow[i];
     10507    if (value<rowLower_[i]-allowedError||value>rowUpper_[i]+allowedError)
     10508      nBad++;
     10509  }
     10510  if (!nBad) {
     10511    // replace
     10512    memcpy(columnLower_,tempColumn,numberColumns_*sizeof(double));
     10513    memcpy(rowLower_,tempRow,numberRows_*sizeof(double));
     10514  }
     10515  delete [] tempColumn;
     10516  return nBad;
     10517}
    1046610518#ifndef SLIM_CLP
    1046710519#include "CoinWarmStartBasis.hpp"
  • trunk/Clp/src/ClpSimplex.hpp

    r2141 r2146  
    360360     */
    361361     int cleanup(int cleanupScaling);
     362     /** Clean primal solution
     363         If you expect solution to only have exact multiples of "exactMultiple" then
     364         this tries moving solution values to nearest multiple.  If still feasible
     365         then the solution is replaced.
     366
     367         This is designed for the case where values should be integral, but Clp may
     368         have values at e.g. 1.0e-13
     369         Returns 0 if successful, n if n rhs violated
     370         The dual version may be written if this gets used.
     371      */
     372      int cleanPrimalSolution(double exactMultiple);
    362373     /** Dual ranging.
    363374         This computes increase/decrease in cost for each given variable and corresponding
  • trunk/Clp/src/ClpSimplexDual.cpp

    r2140 r2146  
    489489                                ifValuesPass);
    490490          if ((specialOptions_&2097152)!=0&&problemStatus_==1&&!ray_&&
    491               !numberRayTries && numberIterations_) {
     491              !numberRayTries && numberIterations_) {
    492492            numberRayTries=1;
    493493            problemStatus_=-1;
     
    22782278     return returnCode;
    22792279}
     2280//#define ABOCA_LITE 4
     2281#if ABOCA_LITE
     2282#include <cilk/cilk.h>
     2283typedef struct {
     2284  double tolerance;
     2285  double theta;
     2286  double * reducedCost;
     2287  const double * lower;
     2288  const double * upper;
     2289  double * work;
     2290  const unsigned char * statusArray;
     2291  int * which;
     2292  int numberInfeasibilities;
     2293  int numberToDo;
     2294} update_duals;
     2295static void
     2296updateDualBit(update_duals & info)
     2297{
     2298  int numberInfeasibilities = 0;
     2299  double tolerance = info.tolerance;
     2300  double theta = info.theta;
     2301  double * COIN_RESTRICT reducedCost = info.reducedCost;
     2302  const double * COIN_RESTRICT lower = info.lower;
     2303  const double * COIN_RESTRICT upper = info.upper;
     2304  double * COIN_RESTRICT work = info.work;
     2305  int number = info.numberToDo;
     2306  int * COIN_RESTRICT which = info.which;
     2307  const unsigned char * COIN_RESTRICT statusArray = info.statusArray;
     2308  double multiplier[] = { -1.0, 1.0};
     2309  for (int i = 0; i < number; i++) {
     2310    int iSequence = which[i];
     2311    double alphaI = work[i];
     2312    work[i] = 0.0;
     2313   
     2314    int iStatus = (statusArray[iSequence] & 3) - 1;
     2315    if (iStatus) {
     2316      double value = reducedCost[iSequence] - theta * alphaI;
     2317      reducedCost[iSequence] = value;
     2318      //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
     2319      double mult = multiplier[iStatus-1];
     2320      value *= mult;
     2321      // skip if free
     2322      if (value < -tolerance&&iStatus > 0) {
     2323        // flipping bounds
     2324        double movement = mult * (upper[iSequence] - lower[iSequence]);
     2325        work[numberInfeasibilities] = movement;
     2326        which[numberInfeasibilities++] = iSequence;
     2327      }
     2328    }
     2329  }
     2330  info.numberInfeasibilities=numberInfeasibilities;
     2331}
     2332#endif
    22802333/* The duals are updated by the given arrays.
    22812334   Returns number of infeasibilities.
     
    23332386                    double mult = multiplier[iStatus-1];
    23342387                    value *= mult;
    2335                     // skip free
    2336                     if (value < -tolerance && iStatus > 0) {
     2388                    // skip if free
     2389                    if (value < -tolerance&&iStatus > 0) {
    23372390                         // flipping bounds
    23382391                         double movement = mult * (lower[iSequence] - upper[iSequence]);
     
    23662419          if ((moreSpecialOptions_ & 8) != 0) {
    23672420               const unsigned char * COIN_RESTRICT statusArray = status_;
     2421#if ABOCA_LITE
     2422               update_duals info[ABOCA_LITE];
     2423               int chunk = (number+ABOCA_LITE-1)/ABOCA_LITE;
     2424               int n=0;
     2425               int * whichX = which;
     2426               for (i=0;i<ABOCA_LITE;i++) {
     2427                 info[i].theta=theta;
     2428                 info[i].tolerance=tolerance;
     2429                 info[i].reducedCost = reducedCost;
     2430                 info[i].lower = lower;
     2431                 info[i].upper = upper;
     2432                 info[i].statusArray=statusArray;
     2433                 info[i].which=which+n;
     2434                 info[i].work=work+n;
     2435                 info[i].numberToDo=CoinMin(chunk,number-n);
     2436                 n += chunk;
     2437               }
     2438               for (i=0;i<ABOCA_LITE;i++) {
     2439                 cilk_spawn updateDualBit(info[i]);
     2440               }
     2441               cilk_sync;
     2442               for (i=0;i<ABOCA_LITE;i++) {
     2443                 int n = info[i].numberInfeasibilities;
     2444                 double * workV = info[i].work;
     2445                 int * whichV = info[i].which;
     2446                 for (int j = 0; j < n; j++) {
     2447                   int iSequence = whichV[j];
     2448                   double movement = workV[j];
     2449                   workV[j] = 0.0;
     2450                   whichX[numberInfeasibilities++]=iSequence;
     2451#ifndef NDEBUG
     2452                   if (fabs(movement) >= 1.0e30)
     2453                     resetFakeBounds(-1000 - iSequence);
     2454#endif
     2455                    changeObj += movement * cost[iSequence];
     2456                    matrix_->add(this, outputArray, iSequence, movement);
     2457                 }
     2458               }
     2459#else
    23682460               for (i = 0; i < number; i++) {
    23692461                    int iSequence = which[i];
     
    23742466                    if (iStatus) {
    23752467                         double value = reducedCost[iSequence] - theta * alphaI;
    2376                          reducedCost[iSequence] = value;
     2468                         reducedCost[iSequence] = value;
     2469                         //printf("xx %d %.18g\n",iSequence,reducedCost[iSequence]);
    23772470                         double mult = multiplier[iStatus-1];
    23782471                         value *= mult;
    2379                          // skip free
    2380                          if (value < -tolerance && iStatus > 0) {
     2472                         // skip if free
     2473                         if (value < -tolerance&&iStatus > 0) {
    23812474                              // flipping bounds
    23822475                              double movement = mult * (upper[iSequence] - lower[iSequence]);
     
    23962489                    }
    23972490               }
     2491#endif
    23982492          } else {
    23992493               for (i = 0; i < number; i++) {
     
    32363330     }
    32373331}
     3332#if ABOCA_LITE
     3333typedef struct {
     3334  const int * COIN_RESTRICT which;
     3335  const double * COIN_RESTRICT work;
     3336  int * COIN_RESTRICT index;
     3337  double * COIN_RESTRICT spare;
     3338  const unsigned char * COIN_RESTRICT status;
     3339  const double * COIN_RESTRICT reducedCost;
     3340  double upperTheta;
     3341  double bestPossible;
     3342  double acceptablePivot;
     3343  double dualTolerance;
     3344  int numberRemaining;
     3345  int numberToDo;
     3346} pricingInfo;
     3347
     3348/* Meat of transposeTimes by column when not scaled and skipping
     3349   and doing part of dualColumn */
     3350static void
     3351dualColumn00(pricingInfo & info)
     3352{
     3353  const int * COIN_RESTRICT which = info.which;
     3354  const double * COIN_RESTRICT work = info.work;
     3355  int * COIN_RESTRICT index = info.index;
     3356  double * COIN_RESTRICT spare = info.spare;
     3357  const unsigned char * COIN_RESTRICT status = info.status;
     3358  const double * COIN_RESTRICT reducedCost = info.reducedCost;
     3359  double upperTheta = info.upperTheta;
     3360  double acceptablePivot = info.acceptablePivot;
     3361  double dualTolerance = info.dualTolerance;
     3362  double bestPossible = info.bestPossible;
     3363  int numberToDo=info.numberToDo;
     3364  double tentativeTheta = 1.0e15;
     3365  int numberRemaining = 0;
     3366  double multiplier[] = { -1.0, 1.0};
     3367  double dualT = - dualTolerance;
     3368  for (int i = 0; i < numberToDo; i++) {
     3369    int iSequence = which[i];
     3370    int wanted = (status[iSequence] & 3) - 1;
     3371    if (wanted) {
     3372      double mult = multiplier[wanted-1];
     3373      double alpha = work[i] * mult;
     3374      if (alpha > 0.0) {
     3375        double oldValue = reducedCost[iSequence] * mult;
     3376        double value = oldValue - tentativeTheta * alpha;
     3377        if (value < dualT) {
     3378          bestPossible = CoinMax(bestPossible, alpha);
     3379          value = oldValue - upperTheta * alpha;
     3380          if (value < dualT && alpha >= acceptablePivot) {
     3381            upperTheta = (oldValue - dualT) / alpha;
     3382          }
     3383          // add to list
     3384          spare[numberRemaining] = alpha * mult;
     3385          index[numberRemaining++] = iSequence;
     3386        }
     3387      }
     3388    }
     3389  }
     3390  info.numberRemaining = numberRemaining;
     3391  info.upperTheta = upperTheta;
     3392  info.bestPossible = bestPossible;
     3393}
     3394// later do so less zeroing in first blocks
     3395// and some of it combined for loop to move and zero
     3396static void moveAndZero(double * to, double * from, int n)
     3397{
     3398  long int distance = from-to;
     3399  assert (distance>=0);
     3400  if (distance==0)
     3401    return;
     3402  memmove(to,from,n*sizeof(double));
     3403  if (n<distance) {
     3404    // no overlap
     3405    memset(from,0,n*sizeof(double));
     3406  } else {
     3407    //memmove(to,from,n*sizeof(double));
     3408    memset(to+n,0,distance*sizeof(double));
     3409  }
     3410}
     3411#endif
    32383412int
    32393413ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
     
    32623436          double multiplier[] = { -1.0, 1.0};
    32633437          double dualT = - dualTolerance_;
    3264           for (int iSection = 0; iSection < 2; iSection++) {
     3438#if ABOCA_LITE==0
     3439          int nSections=2;
     3440#else
     3441          int nSections=1;
     3442#endif
     3443          for (int iSection = 0; iSection < nSections; iSection++) {
    32653444
    32663445               int addSequence;
     
    33123491               }
    33133492          }
     3493#if ABOCA_LITE
     3494          work = columnArray->denseVector();
     3495          number = columnArray->getNumElements();
     3496          which = columnArray->getIndices();
     3497          reducedCost = reducedCostWork_;
     3498          unsigned char * statusArray = status_;
     3499         
     3500          pricingInfo info[ABOCA_LITE];
     3501          int chunk = (number+ABOCA_LITE-1)/ABOCA_LITE;
     3502          int n=0;
     3503          int nR=numberRemaining;
     3504          for (int i=0;i<ABOCA_LITE;i++) {
     3505            info[i].which=which+n;
     3506            info[i].work=work+n;
     3507            info[i].numberToDo=CoinMin(chunk,number-n);
     3508            n += chunk;
     3509            info[i].index = index+nR;
     3510            info[i].spare = spare+nR;
     3511            nR += chunk;
     3512            info[i].reducedCost = reducedCost;
     3513            info[i].upperTheta = upperTheta;
     3514            info[i].bestPossible = bestPossible;
     3515            info[i].acceptablePivot = acceptablePivot;
     3516            info[i].status = statusArray;
     3517            info[i].dualTolerance=dualTolerance_;
     3518          }
     3519          for (int i=0;i<ABOCA_LITE;i++) {
     3520            cilk_spawn dualColumn00(info[i]);
     3521          }
     3522          cilk_sync;
     3523          numberRemaining += info[0].numberRemaining;
     3524          bestPossible = CoinMax(bestPossible,info[0].bestPossible);
     3525          upperTheta = CoinMin(upperTheta,info[0].upperTheta);
     3526          for (int i=1;i<ABOCA_LITE;i++) {
     3527            memmove(index+numberRemaining,info[i].index,info[i].numberRemaining*sizeof(int));
     3528            moveAndZero(spare+numberRemaining,info[i].spare,info[i].numberRemaining);
     3529            numberRemaining += info[i].numberRemaining;
     3530            bestPossible = CoinMax(bestPossible,info[i].bestPossible);
     3531            upperTheta = CoinMin(upperTheta,info[i].upperTheta);
     3532          }
     3533#endif
    33143534     } else {
    33153535          // some free or super basic
  • trunk/Clp/src/ClpSolver.cpp

    r2078 r2146  
    911911                                   if (preSolveFile)
    912912                                        presolveOptions |= 0x40000000;
     913                                   // allow dependency
     914                                   presolveOptions |= 32768;
    913915                                   solveOptions.setPresolveActions(presolveOptions);
    914916                                   solveOptions.setSubstitution(substitution);
  • trunk/Clp/src/CoinAbcBaseFactorization1.cpp

    r2078 r2146  
    24512451#else
    24522452        status=2;
    2453 #if 1 //0 //ABC_NORMAL_DEBUG>1
     2453#if ABC_NORMAL_DEBUG>1
    24542454        std::cout<<"      Went dense at "<<numberRowsLeft_<<" rows "<<
    24552455          totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
Note: See TracChangeset for help on using the changeset viewer.