Changeset 2151


Ignore:
Timestamp:
Jul 9, 2015 5:23:15 AM (3 years ago)
Author:
forrest
Message:

abboca_lite

Location:
stable/1.16/Clp/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • stable/1.16/Clp/src/ClpDualRowSteepest.cpp

    r2070 r2151  
    1111#include "CoinHelperFunctions.hpp"
    1212#include <cstdio>
     13#if defined(ABOCA_LITE) || defined(CILK_OVERLAP)
     14#include <cilk/cilk.h>
     15#endif
    1316//#############################################################################
    1417// Constructors / Destructor / Assignment
     
    661664}
    662665
     666#if ABOCA_LITE
     667typedef struct {
     668  double tolerance;
     669  double primalRatio;
     670  double changeObj;
     671  const double * cost;
     672  double * solution;
     673  const double * lower;
     674  const double * upper;
     675  double * work;
     676  int * which;
     677  double * infeas;
     678  const int * pivotVariable;
     679  int numberAdded;
     680  int numberToDo;
     681#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     682  int numberColumns;
     683#endif
     684} update_primals;
     685static void
     686updatePrimalBit(update_primals & info)
     687{
     688  int numberAdded = 0;
     689  double tolerance = info.tolerance;
     690  double primalRatio = info.primalRatio;
     691  double changeObj=0.0;
     692  const double * COIN_RESTRICT costModel = info.cost;
     693  const double * COIN_RESTRICT lowerModel = info.lower;
     694  const double * COIN_RESTRICT upperModel = info.upper;
     695  double * COIN_RESTRICT solution = info.solution;
     696  double * COIN_RESTRICT infeas = info.infeas;
     697  const int * COIN_RESTRICT pivotVariable = info.pivotVariable;
     698  double * COIN_RESTRICT work = info.work;
     699  int number = info.numberToDo;
     700  int * COIN_RESTRICT which = info.which;
     701#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     702  int numberColumns = info.numberColumns;
     703#endif
     704  for (int i = 0; i < number; i++) {
     705    int iRow = which[i];
     706    int iPivot = pivotVariable[iRow];
     707    double value = solution[iPivot];
     708    double cost = costModel[iPivot];
     709    double change = primalRatio * work[i];
     710    work[i] = 0.0;
     711    value -= change;
     712    changeObj -= change * cost;
     713    double lower = lowerModel[iPivot];
     714    double upper = upperModel[iPivot];
     715    solution[iPivot] = value;
     716    if (value < lower - tolerance) {
     717      value -= lower;
     718      value *= value;
     719#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     720      if (iPivot < numberColumns)
     721        value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     722#endif
     723#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
     724      if (lower == upper)
     725        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     726#endif
     727      // store square in list
     728      if (!infeas[iRow])
     729        which[numberAdded++]=iRow;
     730      infeas[iRow] = value;
     731    } else if (value > upper + tolerance) {
     732      value -= upper;
     733      value *= value;
     734#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     735      if (iPivot < numberColumns)
     736        value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
     737#endif
     738#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
     739      if (lower == upper)
     740        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
     741#endif
     742      // store square in list
     743      if (!infeas[iRow])
     744        which[numberAdded++]=iRow;
     745      infeas[iRow] = value;
     746    } else {
     747      // feasible - was it infeasible - if so set tiny
     748      if (infeas[iRow])
     749        infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
     750    }
     751  }
     752  info.numberAdded = numberAdded;
     753  info.changeObj = changeObj;
     754}
     755#endif
    663756/* Updates primal solution (and maybe list of candidates)
    664757   Uses input vector which it deletes
     
    687780#endif
    688781     if (primalUpdate->packedMode()) {
     782#if ABOCA_LITE == 0
    689783          for (i = 0; i < number; i++) {
    690784               int iRow = which[i];
     
    737831               }
    738832          }
     833#else
     834               update_primals info[ABOCA_LITE];
     835               int chunk = (number+ABOCA_LITE-1)/ABOCA_LITE;
     836               int n=0;
     837               for (i=0;i<ABOCA_LITE;i++) {
     838                 info[i].primalRatio=primalRatio;
     839                 info[i].tolerance=tolerance;
     840                 info[i].cost = costModel;
     841                 info[i].solution = solution;
     842                 info[i].lower = lowerModel;
     843                 info[i].upper = upperModel;
     844                 info[i].infeas = infeas;
     845                 info[i].pivotVariable = pivotVariable;
     846                 info[i].which=which+n;
     847                 info[i].work=work+n;
     848                 info[i].numberToDo=CoinMin(chunk,number-n);
     849#ifdef CLP_DUAL_COLUMN_MULTIPLIER
     850                 info.numberColumns = numberColumns;
     851#endif
     852                 n += chunk;
     853               }
     854               for (i=0;i<ABOCA_LITE;i++) {
     855                 cilk_spawn updatePrimalBit(info[i]);
     856                 //updatePrimalBit(info[i]);
     857               }
     858               cilk_sync;
     859               n=infeasible_->getNumElements();
     860               int * COIN_RESTRICT index = infeasible_->getIndices();
     861               for (i=0;i<ABOCA_LITE;i++) {
     862                 int numberAdded=info[i].numberAdded;
     863                 int * which=info[i].which;
     864                 for (int j=0;j<numberAdded;j++)
     865                   index[n++]=which[j];
     866                 changeObj += info[i].changeObj;
     867               }
     868               infeasible_->setNumElements(n);
     869#endif
    739870     } else {
    740871          for (i = 0; i < number; i++) {
  • stable/1.16/Clp/src/ClpPrimalColumnSteepest.cpp

    r2074 r2151  
    28772877       int sequenceIn = model_->sequenceIn();
    28782878       assert (sequenceIn>=0&&sequenceIn<model_->numberRows()+model_->numberColumns());
    2879        if (weights_[sequenceIn]==(mode_!=1) ? 1.0 : 1.0+ADD_ONE)
     2879       if (weights_[sequenceIn]==((mode_!=1) ? 1.0 : 1.0+ADD_ONE))
    28802880         return;
    28812881       else
  • stable/1.16/Clp/src/ClpSimplexDual.cpp

    r2139 r2151  
    489489                                ifValuesPass);
    490490          if ((specialOptions_&2097152)!=0&&problemStatus_==1&&!ray_&&
    491               !numberRayTries) {
     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.
     
    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;
     
    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
Note: See TracChangeset for help on using the changeset viewer.