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

stuff from stable and start aboca_lite

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.