Changeset 740


Ignore:
Timestamp:
Mar 11, 2006 4:45:45 PM (14 years ago)
Author:
forrest
Message:

makes some problems faster

Location:
trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyWssmp.cpp

    r738 r740  
    299299  int i1=1;
    300300#if WSSMP_BARRIER>=2
    301   wsetmaxthrds(&i1);
     301  int i2=1;
     302  if (model->numberThreads()<=0)
     303    i2=1;
     304  else
     305    i2=model->numberThreads();
     306  wsetmaxthrds(&i2);
    302307#endif
    303308  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
  • trunk/ClpCholeskyWssmpKKT.cpp

    r738 r740  
    234234  int i1=1;
    235235#if WSSMP_BARRIER>=2
    236   wsetmaxthrds(&i1);
     236  int i2=1;
     237  if (model->numberThreads()<=0)
     238    i2=1;
     239  else
     240    i2=model->numberThreads();
     241  wsetmaxthrds(&i2);
    237242#endif
    238243  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
  • trunk/ClpModel.cpp

    r733 r740  
    7676  secondaryStatus_(0),
    7777  lengthNames_(0),
     78  numberThreads_(0),
    7879#ifndef CLP_NO_STD
    7980  defaultHandler_(true),
     
    680681    }
    681682#endif
     683    numberThreads_ = rhs.numberThreads_;
    682684    integerType_ = CoinCopyOfArray(rhs.integerType_,numberColumns_);
    683685    if (rhs.rowActivity_) {
     
    743745    //columnScale_ = rhs.columnScale_;
    744746    lengthNames_ = 0;
     747    numberThreads_ = rhs.numberThreads_;
    745748#ifndef CLP_NO_STD
    746749    rowNames_ = std::vector<std::string> ();
     
    27812784  numberColumns_ = numberColumns;
    27822785  userPointer_ = rhs->userPointer_;
     2786  numberThreads_=0;
    27832787#ifndef CLP_NO_STD
    27842788  if (!dropNames) {
  • trunk/ClpPackedMatrix.cpp

    r683 r740  
    1010
    1111#include "ClpSimplex.hpp"
     12#include "ClpSimplexDual.hpp"
    1213#include "ClpFactorization.hpp"
    1314#ifndef SLIM_CLP
    1415#include "ClpQuadraticObjective.hpp"
    1516#endif
     17//#define THREAD
    1618// at end to get min/max!
    1719#include "ClpPackedMatrix.hpp"
    1820#include "ClpMessage.hpp"
     21//#define COIN_PREFETCH
     22#ifdef COIN_PREFETCH
     23#if 1
     24#define coin_prefetch(mem) \
     25         __asm__ __volatile__ ("prefetchnta %0" : : "m" (*((char *)(mem))))
     26#else
     27#define coin_prefetch(mem) \
     28         __asm__ __volatile__ ("prefetch %0" : : "m" (*((char *)(mem))))
     29#endif
     30#else
     31// dummy
     32#define coin_prefetch(mem)
     33#endif
    1934
    2035//#############################################################################
     
    3045    numberActiveColumns_(0),
    3146    zeroElements_(false),
    32     hasGaps_(true)
     47    hasGaps_(true),
     48    rowCopy_(NULL)
    3349{
    3450  setType(1);
     
    5167    rhsOffset_=NULL;
    5268  }
    53  
     69  if (rhs.rowCopy_) {
     70    rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
     71  } else {
     72    rowCopy_ = NULL;
     73  }
    5474}
    5575
     
    6484  hasGaps_ = true;
    6585  numberActiveColumns_ = matrix_->getNumCols();
     86  rowCopy_ = NULL;
    6687  setType(1);
    6788 
     
    7394  matrix_ = new CoinPackedMatrix(rhs);
    7495  numberActiveColumns_ = matrix_->getNumCols();
     96  rowCopy_ = NULL;
    7597  zeroElements_ = false;
    7698  hasGaps_ = true;
     
    85107{
    86108  delete matrix_;
     109  delete rowCopy_;
    87110}
    88111
     
    100123    zeroElements_ = rhs.zeroElements_;
    101124    hasGaps_ = rhs.hasGaps_;
     125    delete rowCopy_;
     126    if (rhs.rowCopy_) {
     127      rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
     128    } else {
     129      rowCopy_ = NULL;
     130    }
    102131  }
    103132  return *this;
     
    133162  zeroElements_ = rhs.zeroElements_;
    134163  hasGaps_ = rhs.hasGaps_;
     164  rowCopy_ = NULL;
    135165}
    136166ClpPackedMatrix::ClpPackedMatrix (
     
    145175  zeroElements_ = false;
    146176  hasGaps_=true;
     177  rowCopy_ = NULL;
    147178  setType(1);
    148179}
     
    406437    factor *= 0.9;
    407438  assert (!y->getNumElements());
     439  double multiplierX=0.8;
     440  double factor2 = factor*multiplierX;
     441  if (packed&&rowCopy_&&numberInRowArray>2&&numberInRowArray>factor2*numberRows&&
     442      numberInRowArray<0.9*numberRows) {
     443    rowCopy_->transposeTimes(model,rowCopy->getPackedMatrix(),rowArray,y,columnArray);
     444    return;
     445  }
    408446  if (numberInRowArray>factor*numberRows||!rowCopy) {
    409447    // do by column
     
    28282866  return numberErrors;
    28292867}
     2868void
     2869ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
     2870{
     2871  delete rowCopy_;
     2872  rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
     2873  // See if anything in it
     2874  if (!rowCopy_->usefulInfo()) {
     2875    delete rowCopy_;
     2876    rowCopy_=NULL;
     2877  }
     2878}
     2879//#############################################################################
     2880// Constructors / Destructor / Assignment
     2881//#############################################################################
     2882
     2883//-------------------------------------------------------------------
     2884// Default Constructor
     2885//-------------------------------------------------------------------
     2886ClpPackedMatrix2::ClpPackedMatrix2 ()
     2887  : numberBlocks_(0),
     2888    numberRows_(0),
     2889    offset_(NULL),
     2890    count_(NULL),
     2891    rowStart_(NULL),
     2892    column_(NULL),
     2893    work_(NULL)
     2894{
     2895#ifdef THREAD
     2896  threadId_ = NULL;
     2897  info_ = NULL;
     2898#endif
     2899}
     2900//-------------------------------------------------------------------
     2901// Useful Constructor
     2902//-------------------------------------------------------------------
     2903ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy)
     2904  : numberBlocks_(0),
     2905    numberRows_(0),
     2906    offset_(NULL),
     2907    count_(NULL),
     2908    rowStart_(NULL),
     2909    column_(NULL),
     2910    work_(NULL)
     2911{
     2912#ifdef THREAD
     2913  threadId_ = NULL;
     2914  info_ = NULL;
     2915#endif
     2916  numberRows_ = rowCopy->getNumRows();
     2917  if (!numberRows_)
     2918    return;
     2919  int numberColumns = rowCopy->getNumCols();
     2920  const int * column = rowCopy->getIndices();
     2921  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     2922  const int * length = rowCopy->getVectorLengths();
     2923  const double * element = rowCopy->getElements();
     2924  int chunk = 32768; // tune
     2925  //chunk=100;
     2926  // tune
     2927#if 0
     2928  int chunkY[5]={4096,8192,16384,32768,65535};
     2929  int its = model->maximumIterations();
     2930  if (its>=1000000&&its<1000999) {
     2931    its -= 1000000;
     2932    its = its/10;
     2933    if (its>=5) abort();
     2934    chunk=chunkY[its];
     2935    printf("chunk size %d\n",chunk);
     2936#define cpuid(func,ax,bx,cx,dx)\
     2937    __asm__ __volatile__ ("cpuid":\
     2938    "=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
     2939    unsigned int a,b,c,d;
     2940    int func=0;
     2941    cpuid(func,a,b,c,d);
     2942    {
     2943      int i;
     2944      unsigned int value;
     2945      value=b;
     2946      for (i=0;i<4;i++) {
     2947        printf("%c",(value&0xff));
     2948        value = value >>8;
     2949      }
     2950      value=d;
     2951      for (i=0;i<4;i++) {
     2952        printf("%c",(value&0xff));
     2953        value = value >>8;
     2954      }
     2955      value=c;
     2956      for (i=0;i<4;i++) {
     2957        printf("%c",(value&0xff));
     2958        value = value >>8;
     2959      }
     2960      printf("\n");
     2961      int maxfunc = a;
     2962      if (maxfunc>10) {
     2963        printf("not intel?\n");
     2964        abort();
     2965      }
     2966      for (func=1;func<=maxfunc;func++) {
     2967        cpuid(func,a,b,c,d);
     2968        printf("func %d, %x %x %x %x\n",func,a,b,c,d);
     2969      }
     2970    }
     2971#else
     2972    if (numberColumns>10000||chunk==100) {
     2973#endif
     2974  } else {
     2975    //printf("no chunk\n");
     2976    return;
     2977  }
     2978  // Could also analyze matrix to get natural breaks
     2979  numberBlocks_ = (numberColumns+chunk-1)/chunk;
     2980#ifdef THREAD
     2981  // Get work areas
     2982  threadId_ = new pthread_t [numberBlocks_];
     2983  info_ = new dualColumn0Struct[numberBlocks_];
     2984#endif
     2985  // Even out
     2986  chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
     2987  offset_ = new int[numberBlocks_+1];
     2988  offset_[numberBlocks_]=numberColumns;
     2989  int nRow = numberBlocks_*numberRows_;
     2990  count_ = new unsigned short[nRow];
     2991  memset(count_,0,nRow*sizeof(unsigned short));
     2992  rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
     2993  CoinBigIndex nElement = rowStart[numberRows_];
     2994  rowStart_[nRow+numberRows_]=nElement;
     2995  column_ = new unsigned short [nElement];
     2996  // assumes int <= double
     2997  int sizeWork = 6*numberBlocks_;
     2998  work_ = new double[sizeWork];;
     2999  int iBlock;
     3000  int nZero=0;
     3001  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
     3002    int start = iBlock*chunk;
     3003    offset_[iBlock]=start;
     3004    int end = start+chunk;
     3005    for (int iRow=0;iRow<numberRows_;iRow++) {
     3006      if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
     3007        printf("not packed correctly - gaps\n");
     3008        abort();
     3009      }
     3010      bool lastFound=false;
     3011      int nFound=0;
     3012      for (CoinBigIndex j = rowStart[iRow];
     3013           j<rowStart[iRow]+length[iRow];j++) {
     3014        int iColumn = column[j];
     3015        if (iColumn>=start) {
     3016          if (iColumn<end) {
     3017            if (!element[j]) {
     3018              printf("not packed correctly - zero element\n");
     3019              abort();
     3020            }
     3021            column_[j]=iColumn-start;
     3022            nFound++;
     3023            if (lastFound) {
     3024              printf("not packed correctly - out of order\n");
     3025              abort();
     3026            }
     3027          } else {
     3028            //can't find any more
     3029            lastFound=true;
     3030          }
     3031        }
     3032      }
     3033      count_[iRow*numberBlocks_+iBlock]=nFound;
     3034      if (!nFound)
     3035        nZero++;
     3036    }
     3037  }
     3038  //double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
     3039  //printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
     3040}
     3041
     3042//-------------------------------------------------------------------
     3043// Copy constructor
     3044//-------------------------------------------------------------------
     3045ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs)
     3046  : numberBlocks_(rhs.numberBlocks_),
     3047    numberRows_(rhs.numberRows_)
     3048
     3049  if (numberBlocks_) {
     3050    offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
     3051    int nRow = numberBlocks_*numberRows_;
     3052    count_ = CoinCopyOfArray(rhs.count_,nRow);
     3053    rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
     3054    CoinBigIndex nElement = rowStart_[nRow+numberRows_];
     3055    column_ = CoinCopyOfArray(rhs.column_,nElement);
     3056    int sizeWork = 6*numberBlocks_;
     3057    work_ = CoinCopyOfArray(rhs.work_,sizeWork);
     3058#ifdef THREAD
     3059    threadId_ = new pthread_t [numberBlocks_];
     3060    info_ = new dualColumn0Struct[numberBlocks_];
     3061#endif
     3062  } else {
     3063    offset_ = NULL;
     3064    count_ = NULL;
     3065    rowStart_ = NULL;
     3066    column_ = NULL;
     3067    work_ = NULL;
     3068#ifdef THREAD
     3069    threadId_ = NULL;
     3070    info_ = NULL;
     3071#endif
     3072  }
     3073}
     3074//-------------------------------------------------------------------
     3075// Destructor
     3076//-------------------------------------------------------------------
     3077ClpPackedMatrix2::~ClpPackedMatrix2 ()
     3078{
     3079  delete [] offset_;
     3080  delete [] count_;
     3081  delete [] rowStart_;
     3082  delete [] column_;
     3083  delete [] work_;
     3084#ifdef THREAD
     3085  delete [] threadId_;
     3086  delete [] info_;
     3087#endif
     3088}
     3089
     3090//----------------------------------------------------------------
     3091// Assignment operator
     3092//-------------------------------------------------------------------
     3093ClpPackedMatrix2 &
     3094ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
     3095{
     3096  if (this != &rhs) {
     3097    numberBlocks_ = rhs.numberBlocks_;
     3098    numberRows_ = rhs.numberRows_;
     3099    delete [] offset_;
     3100    delete [] count_;
     3101    delete [] rowStart_;
     3102    delete [] column_;
     3103    delete [] work_;
     3104#ifdef THREAD
     3105    delete [] threadId_;
     3106    delete [] info_;
     3107#endif
     3108    if (numberBlocks_) {
     3109      offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
     3110      int nRow = numberBlocks_*numberRows_;
     3111      count_ = CoinCopyOfArray(rhs.count_,nRow);
     3112      rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
     3113      CoinBigIndex nElement = rowStart_[nRow+numberRows_];
     3114      column_ = CoinCopyOfArray(rhs.column_,nElement);
     3115      int sizeWork = 6*numberBlocks_;
     3116      work_ = CoinCopyOfArray(rhs.work_,sizeWork);
     3117#ifdef THREAD
     3118      threadId_ = new pthread_t [numberBlocks_];
     3119      info_ = new dualColumn0Struct[numberBlocks_];
     3120#endif
     3121    } else {
     3122      offset_ = NULL;
     3123      count_ = NULL;
     3124      rowStart_ = NULL;
     3125      column_ = NULL;
     3126      work_ = NULL;
     3127#ifdef THREAD
     3128      threadId_ = NULL;
     3129      info_ = NULL;
     3130#endif
     3131    }
     3132  }
     3133  return *this;
     3134}
     3135static int dualColumn0(const ClpSimplex * model,double * spare,
     3136                       int * spareIndex,const double * arrayTemp,
     3137                       const int * indexTemp,int numberIn,
     3138                       int offset, double acceptablePivot,double * bestPossiblePtr,
     3139                       double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
     3140{
     3141  // do dualColumn0
     3142  int i;
     3143  int numberRemaining=0;
     3144  double bestPossible=0.0;
     3145  double upperTheta=1.0e31;
     3146  double freePivot = acceptablePivot;
     3147  int posFree=-1;
     3148  const double * reducedCost=model->djRegion(1);
     3149  double dualTolerance = model->dualTolerance();
     3150  // We can also see if infeasible or pivoting on free
     3151  double tentativeTheta = 1.0e25;
     3152  for (i=0;i<numberIn;i++) {
     3153    double alpha = arrayTemp[i];
     3154    int iSequence = indexTemp[i]+offset;
     3155    double oldValue;
     3156    double value;
     3157    bool keep;
     3158   
     3159    switch(model->getStatus(iSequence)) {
     3160     
     3161    case ClpSimplex::basic:
     3162    case ClpSimplex::isFixed:
     3163      break;
     3164    case ClpSimplex::isFree:
     3165    case ClpSimplex::superBasic:
     3166      bestPossible = CoinMax(bestPossible,fabs(alpha));
     3167      oldValue = reducedCost[iSequence];
     3168      // If free has to be very large - should come in via dualRow
     3169      if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
     3170        break;
     3171      if (oldValue>dualTolerance) {
     3172        keep = true;
     3173      } else if (oldValue<-dualTolerance) {
     3174        keep = true;
     3175      } else {
     3176        if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5))
     3177          keep = true;
     3178        else
     3179          keep=false;
     3180      }
     3181      if (keep) {
     3182        // free - choose largest
     3183        if (fabs(alpha)>freePivot) {
     3184          freePivot=fabs(alpha);
     3185          posFree = i;
     3186        }
     3187      }
     3188      break;
     3189    case ClpSimplex::atUpperBound:
     3190      oldValue = reducedCost[iSequence];
     3191      value = oldValue-tentativeTheta*alpha;
     3192      //assert (oldValue<=dualTolerance*1.0001);
     3193      if (value>dualTolerance) {
     3194        bestPossible = CoinMax(bestPossible,-alpha);
     3195        value = oldValue-upperTheta*alpha;
     3196        if (value>dualTolerance && -alpha>=acceptablePivot)
     3197          upperTheta = (oldValue-dualTolerance)/alpha;
     3198        // add to list
     3199        spare[numberRemaining]=alpha;
     3200        spareIndex[numberRemaining++]=iSequence;
     3201      }
     3202      break;
     3203    case ClpSimplex::atLowerBound:
     3204      oldValue = reducedCost[iSequence];
     3205      value = oldValue-tentativeTheta*alpha;
     3206      //assert (oldValue>=-dualTolerance*1.0001);
     3207      if (value<-dualTolerance) {
     3208        bestPossible = CoinMax(bestPossible,alpha);
     3209        value = oldValue-upperTheta*alpha;
     3210        if (value<-dualTolerance && alpha>=acceptablePivot)
     3211          upperTheta = (oldValue+dualTolerance)/alpha;
     3212        // add to list
     3213        spare[numberRemaining]=alpha;
     3214        spareIndex[numberRemaining++]=iSequence;
     3215      }
     3216      break;
     3217    }
     3218  }
     3219  *bestPossiblePtr = bestPossible;
     3220  *upperThetaPtr = upperTheta;
     3221  *freePivotPtr = freePivot;
     3222  *posFreePtr = posFree;
     3223  return numberRemaining;
     3224}
     3225static int doOneBlock(double * array, int * index,
     3226                      const double * pi,const CoinBigIndex * rowStart, const double * element,
     3227                      const unsigned short * column, int numberInRowArray, int numberLook)
     3228{
     3229  int iWhich=0;
     3230  int nextN=0;
     3231  CoinBigIndex nextStart=0;
     3232  double nextPi=0.0;
     3233  for (;iWhich<numberInRowArray;iWhich++) {
     3234    nextStart = rowStart[0];
     3235    nextN = rowStart[numberInRowArray]-nextStart;
     3236    rowStart++;
     3237    if (nextN) {
     3238      nextPi = pi[iWhich];
     3239      break;
     3240    }
     3241  }
     3242  int i;
     3243  while (iWhich<numberInRowArray) {
     3244    double value = nextPi;
     3245    CoinBigIndex j=  nextStart;
     3246    int n = nextN;
     3247    // get next
     3248    iWhich++;
     3249    for (;iWhich<numberInRowArray;iWhich++) {
     3250      nextStart = rowStart[0];
     3251      nextN = rowStart[numberInRowArray]-nextStart;
     3252      rowStart++;
     3253      if (nextN) {
     3254        coin_prefetch(element+nextStart);
     3255        nextPi = pi[iWhich];
     3256        break;
     3257      }
     3258    }
     3259    CoinBigIndex end =j+n;
     3260    //coin_prefetch(element+rowStart_[i+1]);
     3261    //coin_prefetch(column_+rowStart_[i+1]);
     3262    if (n<100) {
     3263      if ((n&1)!=0) {
     3264        unsigned int jColumn = column[j];
     3265        array[jColumn] -= value*element[j];
     3266        j++;
     3267      }     
     3268      coin_prefetch(column+nextStart);
     3269      for (;j<end;j+=2) {
     3270        unsigned int jColumn0 = column[j];
     3271        double value0 = value*element[j];
     3272        unsigned int jColumn1 = column[j+1];
     3273        double value1 = value*element[j+1];
     3274        array[jColumn0] -= value0;
     3275        array[jColumn1] -= value1;
     3276      }
     3277    } else {
     3278      if ((n&1)!=0) {
     3279        unsigned int jColumn = column[j];
     3280        array[jColumn] -= value*element[j];
     3281        j++;
     3282      }     
     3283      if ((n&2)!=0) {
     3284        unsigned int jColumn0 = column[j];
     3285        double value0 = value*element[j];
     3286        unsigned int jColumn1 = column[j+1];
     3287        double value1 = value*element[j+1];
     3288        array[jColumn0] -= value0;
     3289        array[jColumn1] -= value1;
     3290        j+=2;
     3291      }     
     3292      if ((n&4)!=0) {
     3293        unsigned int jColumn0 = column[j];
     3294        double value0 = value*element[j];
     3295        unsigned int jColumn1 = column[j+1];
     3296        double value1 = value*element[j+1];
     3297        unsigned int jColumn2 = column[j+2];
     3298        double value2 = value*element[j+2];
     3299        unsigned int jColumn3 = column[j+3];
     3300        double value3 = value*element[j+3];
     3301        array[jColumn0] -= value0;
     3302        array[jColumn1] -= value1;
     3303        array[jColumn2] -= value2;
     3304        array[jColumn3] -= value3;
     3305        j+=4;
     3306      }
     3307      //coin_prefetch(column+nextStart);
     3308      for (;j<end;j+=8) {
     3309        coin_prefetch(element+j+16);
     3310        unsigned int jColumn0 = column[j];
     3311        double value0 = value*element[j];
     3312        unsigned int jColumn1 = column[j+1];
     3313        double value1 = value*element[j+1];
     3314        unsigned int jColumn2 = column[j+2];
     3315        double value2 = value*element[j+2];
     3316        unsigned int jColumn3 = column[j+3];
     3317        double value3 = value*element[j+3];
     3318        array[jColumn0] -= value0;
     3319        array[jColumn1] -= value1;
     3320        array[jColumn2] -= value2;
     3321        array[jColumn3] -= value3;
     3322        coin_prefetch(column+j+16);
     3323        jColumn0 = column[j+4];
     3324        value0 = value*element[j+4];
     3325        jColumn1 = column[j+5];
     3326        value1 = value*element[j+5];
     3327        jColumn2 = column[j+6];
     3328        value2 = value*element[j+6];
     3329        jColumn3 = column[j+7];
     3330        value3 = value*element[j+7];
     3331        array[jColumn0] -= value0;
     3332        array[jColumn1] -= value1;
     3333        array[jColumn2] -= value2;
     3334        array[jColumn3] -= value3;
     3335      }
     3336    }
     3337  }
     3338  // get rid of tiny values
     3339  int nSmall = numberLook;
     3340  int numberNonZero=0;
     3341  for (i=0;i<nSmall;i++) {
     3342    double value = array[i];
     3343    array[i]=0.0;
     3344    if (fabs(value)>1.0e-12) {
     3345      array[numberNonZero]=value;
     3346      index[numberNonZero++]=i;
     3347    }
     3348  }
     3349  for (;i<numberLook;i+=4) {
     3350    double value0 = array[i+0];
     3351    double value1 = array[i+1];
     3352    double value2 = array[i+2];
     3353    double value3 = array[i+3];
     3354    array[i+0]=0.0;
     3355    array[i+1]=0.0;
     3356    array[i+2]=0.0;
     3357    array[i+3]=0.0;
     3358    if (fabs(value0)>1.0e-12) {
     3359      array[numberNonZero]=value0;
     3360      index[numberNonZero++]=i+0;
     3361    }
     3362    if (fabs(value1)>1.0e-12) {
     3363      array[numberNonZero]=value1;
     3364      index[numberNonZero++]=i+1;
     3365    }
     3366    if (fabs(value2)>1.0e-12) {
     3367      array[numberNonZero]=value2;
     3368      index[numberNonZero++]=i+2;
     3369    }
     3370    if (fabs(value3)>1.0e-12) {
     3371      array[numberNonZero]=value3;
     3372      index[numberNonZero++]=i+3;
     3373    }
     3374  }
     3375  return numberNonZero;
     3376}
     3377#ifdef THREAD
     3378static void * doOneBlockThread(void * voidInfo)
     3379{
     3380  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
     3381  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
     3382                                  info->rowStart,info->element,info->column,
     3383                                  info->numberInRowArray,info->numberLook);
     3384  return NULL;
     3385}
     3386static void * doOneBlockAnd0Thread(void * voidInfo)
     3387{
     3388  dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
     3389  *(info->numberInPtr) =  doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
     3390                                  info->rowStart,info->element,info->column,
     3391                                  info->numberInRowArray,info->numberLook);
     3392  *(info->numberOutPtr) =  dualColumn0(info->model,info->spare,
     3393                                    info->spareIndex,(const double *)info->arrayTemp,
     3394                                    (const int *) info->indexTemp,*(info->numberInPtr),
     3395                                    info->offset, info->acceptablePivot,info->bestPossiblePtr,
     3396                                    info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
     3397  return NULL;
     3398}
     3399#endif
     3400/* Return <code>x * scalar * A in <code>z</code>.
     3401   Note - x packed and z will be packed mode
     3402   Squashes small elements and knows about ClpSimplex */
     3403void
     3404ClpPackedMatrix2::transposeTimes(const ClpSimplex * model,
     3405                                 const CoinPackedMatrix * rowCopy,
     3406                                 const CoinIndexedVector * rowArray,
     3407                                 CoinIndexedVector * spareArray,
     3408                                 CoinIndexedVector * columnArray) const
     3409{
     3410  // See if dualColumn0 coding wanted
     3411  bool dualColumn = model->spareIntArray_[0]==1;
     3412  double acceptablePivot = model->spareDoubleArray_[0];
     3413  double bestPossible=0.0;
     3414  double upperTheta=1.0e31;
     3415  double freePivot = acceptablePivot;
     3416  int posFree=-1;
     3417  int numberRemaining=0;
     3418  //if (model->numberIterations()>=200000) {
     3419  //printf("time %g\n",CoinCpuTime()-startTime);
     3420  //exit(77);
     3421  //}
     3422  double * pi = rowArray->denseVector();
     3423  int numberNonZero=0;
     3424  int * index = columnArray->getIndices();
     3425  double * array = columnArray->denseVector();
     3426  int numberInRowArray = rowArray->getNumElements();
     3427  const int * whichRow = rowArray->getIndices();
     3428  double * element = const_cast<double *>(rowCopy->getElements());
     3429  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     3430  int i;
     3431  CoinBigIndex * rowStart2 = rowStart_;
     3432  if (!dualColumn) {
     3433    for (i=0;i<numberInRowArray;i++) {
     3434      int iRow = whichRow[i];
     3435      CoinBigIndex start = rowStart[iRow];
     3436      *rowStart2 = start;
     3437      unsigned short * count1 = count_+iRow*numberBlocks_;
     3438      int put=0;
     3439      for (int j=0;j<numberBlocks_;j++) {
     3440        put += numberInRowArray;
     3441        start += count1[j];
     3442        rowStart2[put]=start;
     3443      }
     3444      rowStart2 ++;
     3445    }
     3446  } else {
     3447    // also do dualColumn stuff
     3448    double * spare = spareArray->denseVector();
     3449    int * spareIndex = spareArray->getIndices();
     3450    const double * reducedCost=model->djRegion(0);
     3451    double dualTolerance = model->dualTolerance();
     3452    // We can also see if infeasible or pivoting on free
     3453    double tentativeTheta = 1.0e25;
     3454    int addSequence=model->numberColumns();
     3455    for (i=0;i<numberInRowArray;i++) {
     3456      int iRow = whichRow[i];
     3457      double alpha=pi[i];
     3458      double oldValue;
     3459      double value;
     3460      bool keep;
     3461
     3462      switch(model->getStatus(iRow+addSequence)) {
     3463         
     3464      case ClpSimplex::basic:
     3465      case ClpSimplex::isFixed:
     3466        break;
     3467      case ClpSimplex::isFree:
     3468      case ClpSimplex::superBasic:
     3469        bestPossible = CoinMax(bestPossible,fabs(alpha));
     3470        oldValue = reducedCost[iRow];
     3471        // If free has to be very large - should come in via dualRow
     3472        if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
     3473          break;
     3474        if (oldValue>dualTolerance) {
     3475          keep = true;
     3476        } else if (oldValue<-dualTolerance) {
     3477          keep = true;
     3478        } else {
     3479          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5))
     3480            keep = true;
     3481          else
     3482            keep=false;
     3483        }
     3484        if (keep) {
     3485          // free - choose largest
     3486          if (fabs(alpha)>freePivot) {
     3487            freePivot=fabs(alpha);
     3488            posFree = i + addSequence;
     3489          }
     3490        }
     3491        break;
     3492      case ClpSimplex::atUpperBound:
     3493        oldValue = reducedCost[iRow];
     3494        value = oldValue-tentativeTheta*alpha;
     3495        //assert (oldValue<=dualTolerance*1.0001);
     3496        if (value>dualTolerance) {
     3497          bestPossible = CoinMax(bestPossible,-alpha);
     3498          value = oldValue-upperTheta*alpha;
     3499          if (value>dualTolerance && -alpha>=acceptablePivot)
     3500            upperTheta = (oldValue-dualTolerance)/alpha;
     3501          // add to list
     3502          spare[numberRemaining]=alpha;
     3503          spareIndex[numberRemaining++]=iRow+addSequence;
     3504        }
     3505        break;
     3506      case ClpSimplex::atLowerBound:
     3507        oldValue = reducedCost[iRow];
     3508        value = oldValue-tentativeTheta*alpha;
     3509        //assert (oldValue>=-dualTolerance*1.0001);
     3510        if (value<-dualTolerance) {
     3511          bestPossible = CoinMax(bestPossible,alpha);
     3512          value = oldValue-upperTheta*alpha;
     3513          if (value<-dualTolerance && alpha>=acceptablePivot)
     3514            upperTheta = (oldValue+dualTolerance)/alpha;
     3515          // add to list
     3516          spare[numberRemaining]=alpha;
     3517          spareIndex[numberRemaining++]=iRow+addSequence;
     3518        }
     3519        break;
     3520      }
     3521      CoinBigIndex start = rowStart[iRow];
     3522      *rowStart2 = start;
     3523      unsigned short * count1 = count_+iRow*numberBlocks_;
     3524      int put=0;
     3525      for (int j=0;j<numberBlocks_;j++) {
     3526        put += numberInRowArray;
     3527        start += count1[j];
     3528        rowStart2[put]=start;
     3529      }
     3530      rowStart2 ++;
     3531    }
     3532  }
     3533  double * spare = spareArray->denseVector();
     3534  int * spareIndex = spareArray->getIndices();
     3535  int saveNumberRemaining=numberRemaining;
     3536  int iBlock;
     3537  for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
     3538    double * dwork = work_+6*iBlock;
     3539    int * iwork = (int *) (dwork+3);
     3540    if (!dualColumn) {
     3541#ifndef THREAD
     3542      int offset = offset_[iBlock];
     3543      int offset3 = offset;
     3544      offset = numberNonZero;
     3545      double * arrayTemp = array+offset;
     3546      int * indexTemp = index+offset;
     3547      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
     3548                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
     3549      int number = iwork[0];
     3550      for (i=0;i<number;i++) {
     3551        //double value = arrayTemp[i];
     3552        //arrayTemp[i]=0.0;
     3553        //array[numberNonZero]=value;
     3554        index[numberNonZero++]=indexTemp[i]+offset3;
     3555      }
     3556#else
     3557      int offset = offset_[iBlock];
     3558      double * arrayTemp = array+offset;
     3559      int * indexTemp = index+offset;
     3560      dualColumn0Struct * infoPtr = info_+iBlock;
     3561      infoPtr->arrayTemp=arrayTemp;
     3562      infoPtr->indexTemp=indexTemp;
     3563      infoPtr->numberInPtr=&iwork[0];
     3564      infoPtr->pi=pi;
     3565      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
     3566      infoPtr->element=element;
     3567      infoPtr->column=column_;
     3568      infoPtr->numberInRowArray=numberInRowArray;
     3569      infoPtr->numberLook=offset_[iBlock+1]-offset;
     3570      pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
     3571#endif
     3572    } else {
     3573#ifndef THREAD
     3574      int offset = offset_[iBlock];
     3575      // allow for already saved
     3576      int offset2 = offset + saveNumberRemaining;
     3577      int offset3 = offset;
     3578      offset = numberNonZero;
     3579      offset2 = numberRemaining;
     3580      double * arrayTemp = array+offset;
     3581      int * indexTemp = index+offset;
     3582      iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
     3583                            element,column_,numberInRowArray,offset_[iBlock+1]-offset);
     3584      iwork[1] = dualColumn0(model,spare+offset2,
     3585                                   spareIndex+offset2,
     3586                                   arrayTemp,indexTemp,
     3587                                   iwork[0],offset3,acceptablePivot,
     3588                                   &dwork[0],&dwork[1],&iwork[2],
     3589                                   &dwork[2]);
     3590      int number = iwork[0];
     3591      int numberLook = iwork[1];
     3592#if 1
     3593      numberRemaining += numberLook;
     3594#else
     3595      double * spareTemp = spare+offset2;
     3596      const int * spareIndexTemp = spareIndex+offset2;
     3597      for (i=0;i<numberLook;i++) {
     3598        double value = spareTemp[i];
     3599        spareTemp[i]=0.0;
     3600        spare[numberRemaining] = value;
     3601        spareIndex[numberRemaining++] = spareIndexTemp[i];
     3602      }
     3603#endif
     3604      if (dwork[2]>freePivot) {
     3605        freePivot=dwork[2];
     3606        posFree = iwork[2]+numberNonZero;
     3607      }
     3608      upperTheta =  CoinMin(dwork[1],upperTheta);
     3609      bestPossible = CoinMax(dwork[0],bestPossible);
     3610      for (i=0;i<number;i++) {
     3611        // double value = arrayTemp[i];
     3612        //arrayTemp[i]=0.0;
     3613        //array[numberNonZero]=value;
     3614        index[numberNonZero++]=indexTemp[i]+offset3;
     3615      }
     3616#else
     3617      int offset = offset_[iBlock];
     3618      // allow for already saved
     3619      int offset2 = offset + saveNumberRemaining;
     3620      double * arrayTemp = array+offset;
     3621      int * indexTemp = index+offset;
     3622      dualColumn0Struct * infoPtr = info_+iBlock;
     3623      infoPtr->model=model;
     3624      infoPtr->spare=spare+offset2;
     3625      infoPtr->spareIndex=spareIndex+offset2;
     3626      infoPtr->arrayTemp=arrayTemp;
     3627      infoPtr->indexTemp=indexTemp;
     3628      infoPtr->numberInPtr=&iwork[0];
     3629      infoPtr->offset=offset;
     3630      infoPtr->acceptablePivot=acceptablePivot;
     3631      infoPtr->bestPossiblePtr=&dwork[0];
     3632      infoPtr->upperThetaPtr=&dwork[1];
     3633      infoPtr->posFreePtr=&iwork[2];
     3634      infoPtr->freePivotPtr=&dwork[2];
     3635      infoPtr->numberOutPtr=&iwork[1];
     3636      infoPtr->pi=pi;
     3637      infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
     3638      infoPtr->element=element;
     3639      infoPtr->column=column_;
     3640      infoPtr->numberInRowArray=numberInRowArray;
     3641      infoPtr->numberLook=offset_[iBlock+1]-offset;
     3642      if (iBlock>=2)
     3643        pthread_join(threadId_[iBlock-2],NULL);
     3644      pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
     3645      //pthread_join(threadId_[iBlock],NULL);
     3646#endif
     3647    }
     3648  }
     3649  for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
     3650#ifdef THREAD
     3651    pthread_join(threadId_[iBlock],NULL);
     3652#endif
     3653  }
     3654#ifdef THREAD
     3655  for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
     3656    //pthread_join(threadId_[iBlock],NULL);
     3657    int offset = offset_[iBlock];
     3658    double * dwork = work_+6*iBlock;
     3659    int * iwork = (int *) (dwork+3);
     3660    int number = iwork[0];
     3661    if (dualColumn) {
     3662      // allow for already saved
     3663      int offset2 = offset + saveNumberRemaining;
     3664      int numberLook = iwork[1];
     3665      double * spareTemp = spare+offset2;
     3666      const int * spareIndexTemp = spareIndex+offset2;
     3667      for (i=0;i<numberLook;i++) {
     3668        double value = spareTemp[i];
     3669        spareTemp[i]=0.0;
     3670        spare[numberRemaining] = value;
     3671        spareIndex[numberRemaining++] = spareIndexTemp[i];
     3672      }
     3673      if (dwork[2]>freePivot) {
     3674        freePivot=dwork[2];
     3675        posFree = iwork[2]+numberNonZero;
     3676      }
     3677      upperTheta =  CoinMin(dwork[1],upperTheta);
     3678      bestPossible = CoinMax(dwork[0],bestPossible);
     3679    }
     3680    double * arrayTemp = array+offset;
     3681    const int * indexTemp = index+offset;
     3682    for (i=0;i<number;i++) {
     3683      double value = arrayTemp[i];
     3684      arrayTemp[i]=0.0;
     3685      array[numberNonZero]=value;
     3686      index[numberNonZero++]=indexTemp[i]+offset;
     3687    }
     3688  }
     3689#endif
     3690  columnArray->setNumElements(numberNonZero);
     3691  columnArray->setPackedMode(true);
     3692  if (dualColumn) {
     3693    model->spareDoubleArray_[0]=upperTheta;
     3694    model->spareDoubleArray_[1]=bestPossible;
     3695    // and theta and alpha and sequence
     3696    if (posFree<0) {
     3697      model->spareIntArray_[1]=-1;
     3698    } else {
     3699      const double * reducedCost=model->djRegion(0);
     3700      double alpha;
     3701      int numberColumns=model->numberColumns();
     3702      if (posFree<numberColumns) {
     3703        alpha = columnArray->denseVector()[posFree];
     3704        posFree = columnArray->getIndices()[posFree];
     3705      } else {
     3706        alpha = rowArray->denseVector()[posFree-numberColumns];
     3707        posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
     3708      }
     3709      model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
     3710      model->spareDoubleArray_[3]=alpha;
     3711      model->spareIntArray_[1]=posFree;
     3712    }
     3713    spareArray->setNumElements(numberRemaining);
     3714    // signal done
     3715    model->spareIntArray_[0]=-1;
     3716  }
     3717}
    28303718#ifdef CLP_ALL_ONE_FILE
    28313719#undef reference
  • trunk/ClpSimplex.cpp

    r737 r740  
    123123    columnArray_[i]=NULL;
    124124  }
     125  for (i=0;i<4;i++) {
     126    spareIntArray_[i]=0;
     127    spareDoubleArray_[i]=0.0;
     128  }
    125129  saveStatus_=NULL;
    126130  // get an empty factorization so we can set tolerances etc
     
    229233    rowArray_[i]=NULL;
    230234    columnArray_[i]=NULL;
     235  }
     236  for (i=0;i<4;i++) {
     237    spareIntArray_[i]=0;
     238    spareDoubleArray_[i]=0.0;
    231239  }
    232240  saveStatus_=NULL;
     
    16601668    columnArray_[i]=NULL;
    16611669  }
     1670  for (i=0;i<4;i++) {
     1671    spareIntArray_[i]=0;
     1672    spareDoubleArray_[i]=0.0;
     1673  }
    16621674  saveStatus_=NULL;
    16631675  factorization_ = NULL;
     
    17611773    columnArray_[i]=NULL;
    17621774  }
     1775  for (i=0;i<4;i++) {
     1776    spareIntArray_[i]=0;
     1777    spareDoubleArray_[i]=0.0;
     1778  }
    17631779  saveStatus_=NULL;
    17641780  // get an empty factorization so we can set tolerances etc
     
    19261942  else
    19271943    progress_=NULL;
     1944  for (int i=0;i<4;i++) {
     1945    spareIntArray_[i]=rhs.spareIntArray_[i];
     1946    spareDoubleArray_[i]=rhs.spareDoubleArray_[i];
     1947  }
    19281948  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    19291949  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     
    27132733    } else if (makeRowCopy&&scalingFlag_>0&&!oldMatrix) {
    27142734      matrix_->scaleRowCopy(this);
     2735    }
     2736    // See if we can try for faster row copy
     2737    if (makeRowCopy&&!oldMatrix) {
     2738      ClpPackedMatrix* clpMatrix =
     2739        dynamic_cast< ClpPackedMatrix*>(matrix_);
     2740      if (clpMatrix&&numberThreads_)
     2741        clpMatrix->specialRowCopy(this,rowCopy_);
    27152742    }
    27162743  }
  • trunk/ClpSimplexDual.cpp

    r736 r740  
    10151015        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
    10161016        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     1017        // Allow to do dualColumn0
     1018        if (numberThreads_<-1)
     1019          spareIntArray_[0]=1;
     1020        spareDoubleArray_[0]=acceptablePivot;
     1021        rowArray_[3]->clear();
     1022        sequenceIn_=-1;
    10171023        // put row of tableau in rowArray[0] and columnArray[0]
    10181024        matrix_->transposeTimes(this,-1.0,
    10191025                              rowArray_[0],rowArray_[3],columnArray_[0]);
    10201026        // do ratio test for normal iteration
    1021         bestPossiblePivot = dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
    1022                  rowArray_[3],acceptablePivot,dubiousWeights);
     1027        bestPossiblePivot = dualColumn(rowArray_[0],columnArray_[0],rowArray_[3],
     1028                 columnArray_[1],acceptablePivot,dubiousWeights);
    10231029      } else {
    10241030        // Make sure direction plausible
     
    23882394  }
    23892395}
    2390 /*
    2391    Row array has row part of pivot row (as duals so sign may be switched)
    2392    Column array has column part.
    2393    This chooses pivot column.
    2394    Spare array will be needed when we start getting clever.
    2395    We will check for basic so spare array will never overflow.
    2396    If necessary will modify costs
    2397 */
    2398 double
    2399 ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
    2400                            CoinIndexedVector * columnArray,
     2396int
     2397ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
     2398                           const CoinIndexedVector * columnArray,
    24012399                           CoinIndexedVector * spareArray,
    2402                            CoinIndexedVector * spareArray2,
    24032400                           double acceptablePivot,
    2404                            CoinBigIndex * dubiousWeights)
     2401                           double & upperReturn, double &bestReturn)
    24052402{
    2406   double * work;
     2403  // do first pass to get possibles
     2404  double * spare = spareArray->denseVector();
     2405  int * index = spareArray->getIndices();
     2406  const double * work;
    24072407  int number;
    2408   int * which;
    2409   double * reducedCost;
    2410   int iSection;
    2411 
    2412   sequenceIn_=-1;
    2413   int numberPossiblySwapped=0;
    2414   int numberRemaining=0;
    2415  
    2416   double totalThru=0.0; // for when variables flip
    2417   //double saveAcceptable=acceptablePivot;
    2418   //acceptablePivot=1.0e-9;
    2419 
    2420   double bestEverPivot=acceptablePivot;
    2421   int lastSequence = -1;
    2422   double lastPivot=0.0;
    2423   double upperTheta;
    2424   double newTolerance = dualTolerance_;
    2425   //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
    2426   // will we need to increase tolerance
    2427   bool thisIncrease=false;
    2428   // If we think we need to modify costs (not if something from broad sweep)
    2429   bool modifyCosts=false;
    2430   // Increase in objective due to swapping bounds (may be negative)
    2431   double increaseInObjective=0.0;
    2432 
    2433   // use spareArrays to put ones looked at in
    2434   // we are going to flip flop between
    2435   int iFlip = 0;
    2436   // Possible list of pivots
    2437   int interesting[2];
    2438   // where possible swapped ones are
    2439   int swapped[2];
    2440   // for zeroing out arrays after
    2441   int marker[2][2];
    2442   // pivot elements
    2443   double * array[2], * spare, * spare2;
    2444   // indices
    2445   int * indices[2], * index, * index2;
    2446   spareArray->clear();
    2447   spareArray2->clear();
    2448   array[0] = spareArray->denseVector();
    2449   indices[0] = spareArray->getIndices();
    2450   spare = array[0];
    2451   index = indices[0];
    2452   array[1] = spareArray2->denseVector();
    2453   indices[1] = spareArray2->getIndices();
    2454   int i;
    2455   const double * lower;
    2456   const double * upper;
    2457 
    2458   // initialize lists
    2459   for (i=0;i<2;i++) {
    2460     interesting[i]=0;
    2461     swapped[i]=numberColumns_;
    2462     marker[i][0]=0;
    2463     marker[i][1]=numberColumns_;
    2464   }
    2465 
    2466   /*
    2467     First we get a list of possible pivots.  We can also see if the
    2468     problem looks infeasible or whether we want to pivot in free variable.
    2469     This may make objective go backwards but can only happen a finite
    2470     number of times and I do want free variables basic.
    2471 
    2472     Then we flip back and forth.  At the start of each iteration
    2473     interesting[iFlip] should have possible candidates and swapped[iFlip]
    2474     will have pivots if we decide to take a previous pivot.
    2475     At end of each iteration interesting[1-iFlip] should have
    2476     candidates if we go through this theta and swapped[1-iFlip]
    2477     pivots if we don't go through.
    2478 
    2479     At first we increase theta and see what happens.  We start
    2480     theta at a reasonable guess.  If in right area then we do bit by bit.
    2481 
    2482    */
    2483 
    2484   // do first pass to get possibles
     2408  const int * which;
     2409  const double * reducedCost;
    24852410  // We can also see if infeasible or pivoting on free
    24862411  double tentativeTheta = 1.0e25;
    2487   upperTheta = 1.0e31;
     2412  double upperTheta = 1.0e31;
    24882413  double freePivot = acceptablePivot;
    24892414  double bestPossible=0.0;
    2490 
    2491   for (iSection=0;iSection<2;iSection++) {
     2415  int numberRemaining=0;
     2416  int i;
     2417  for (int iSection=0;iSection<2;iSection++) {
    24922418
    24932419    int addSequence;
    24942420
    24952421    if (!iSection) {
    2496       lower = rowLowerWork_;
    2497       upper = rowUpperWork_;
    24982422      work = rowArray->denseVector();
    24992423      number = rowArray->getNumElements();
     
    25022426      addSequence = numberColumns_;
    25032427    } else {
    2504       lower = columnLowerWork_;
    2505       upper = columnUpperWork_;
    25062428      work = columnArray->denseVector();
    25072429      number = columnArray->getNumElements();
     
    25562478        value = oldValue-tentativeTheta*alpha;
    25572479        //assert (oldValue<=dualTolerance_*1.0001);
    2558         if (value>newTolerance) {
     2480        if (value>dualTolerance_) {
    25592481          bestPossible = CoinMax(bestPossible,-alpha);
    25602482          value = oldValue-upperTheta*alpha;
    2561           if (value>newTolerance && -alpha>=acceptablePivot)
    2562             upperTheta = (oldValue-newTolerance)/alpha;
     2483          if (value>dualTolerance_ && -alpha>=acceptablePivot)
     2484            upperTheta = (oldValue-dualTolerance_)/alpha;
    25632485          // add to list
    25642486          spare[numberRemaining]=alpha;
     
    25712493        value = oldValue-tentativeTheta*alpha;
    25722494        //assert (oldValue>=-dualTolerance_*1.0001);
    2573         if (value<-newTolerance) {
     2495        if (value<-dualTolerance_) {
    25742496          bestPossible = CoinMax(bestPossible,alpha);
    25752497          value = oldValue-upperTheta*alpha;
    2576           if (value<-newTolerance && alpha>=acceptablePivot)
    2577             upperTheta = (oldValue+newTolerance)/alpha;
     2498          if (value<-dualTolerance_ && alpha>=acceptablePivot)
     2499            upperTheta = (oldValue+dualTolerance_)/alpha;
    25782500          // add to list
    25792501          spare[numberRemaining]=alpha;
     
    25842506    }
    25852507  }
     2508  upperReturn = upperTheta;
     2509  bestReturn = bestPossible;
     2510  return numberRemaining;
     2511}
     2512/*
     2513   Row array has row part of pivot row (as duals so sign may be switched)
     2514   Column array has column part.
     2515   This chooses pivot column.
     2516   Spare array will be needed when we start getting clever.
     2517   We will check for basic so spare array will never overflow.
     2518   If necessary will modify costs
     2519*/
     2520double
     2521ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
     2522                           CoinIndexedVector * columnArray,
     2523                           CoinIndexedVector * spareArray,
     2524                           CoinIndexedVector * spareArray2,
     2525                           double acceptablePivot,
     2526                           CoinBigIndex * dubiousWeights)
     2527{
     2528  int numberPossiblySwapped=0;
     2529  int numberRemaining=0;
     2530 
     2531  double totalThru=0.0; // for when variables flip
     2532  //double saveAcceptable=acceptablePivot;
     2533  //acceptablePivot=1.0e-9;
     2534
     2535  double bestEverPivot=acceptablePivot;
     2536  int lastSequence = -1;
     2537  double lastPivot=0.0;
     2538  double upperTheta;
     2539  double newTolerance = dualTolerance_;
     2540  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
     2541  // will we need to increase tolerance
     2542  bool thisIncrease=false;
     2543  // If we think we need to modify costs (not if something from broad sweep)
     2544  bool modifyCosts=false;
     2545  // Increase in objective due to swapping bounds (may be negative)
     2546  double increaseInObjective=0.0;
     2547
     2548  // use spareArrays to put ones looked at in
     2549  // we are going to flip flop between
     2550  int iFlip = 0;
     2551  // Possible list of pivots
     2552  int interesting[2];
     2553  // where possible swapped ones are
     2554  int swapped[2];
     2555  // for zeroing out arrays after
     2556  int marker[2][2];
     2557  // pivot elements
     2558  double * array[2], * spare, * spare2;
     2559  // indices
     2560  int * indices[2], * index, * index2;
     2561  spareArray2->clear();
     2562  array[0] = spareArray->denseVector();
     2563  indices[0] = spareArray->getIndices();
     2564  spare = array[0];
     2565  index = indices[0];
     2566  array[1] = spareArray2->denseVector();
     2567  indices[1] = spareArray2->getIndices();
     2568  int i;
     2569
     2570  // initialize lists
     2571  for (i=0;i<2;i++) {
     2572    interesting[i]=0;
     2573    swapped[i]=numberColumns_;
     2574    marker[i][0]=0;
     2575    marker[i][1]=numberColumns_;
     2576  }
     2577  /*
     2578    First we get a list of possible pivots.  We can also see if the
     2579    problem looks infeasible or whether we want to pivot in free variable.
     2580    This may make objective go backwards but can only happen a finite
     2581    number of times and I do want free variables basic.
     2582
     2583    Then we flip back and forth.  At the start of each iteration
     2584    interesting[iFlip] should have possible candidates and swapped[iFlip]
     2585    will have pivots if we decide to take a previous pivot.
     2586    At end of each iteration interesting[1-iFlip] should have
     2587    candidates if we go through this theta and swapped[1-iFlip]
     2588    pivots if we don't go through.
     2589
     2590    At first we increase theta and see what happens.  We start
     2591    theta at a reasonable guess.  If in right area then we do bit by bit.
     2592
     2593   */
     2594
     2595  // do first pass to get possibles
     2596  upperTheta = 1.0e31;
     2597  double bestPossible=0.0;
     2598  if (spareIntArray_[0]!=-1) {
     2599    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
     2600                                  acceptablePivot,upperTheta,bestPossible);
     2601  } else {
     2602    // already done
     2603    numberRemaining = spareArray->getNumElements();
     2604    spareArray->setNumElements(0);
     2605    upperTheta = spareDoubleArray_[0];
     2606    bestPossible = spareDoubleArray_[1];
     2607    theta_ = spareDoubleArray_[2];
     2608    alpha_ = spareDoubleArray_[3];
     2609    sequenceIn_ = spareIntArray_[1];
     2610  }
     2611  // switch off
     2612  spareIntArray_[0]=0;
     2613  // We can also see if infeasible or pivoting on free
     2614  double tentativeTheta = 1.0e25;
    25862615  interesting[0]=numberRemaining;
    25872616  marker[0][0] = numberRemaining;
     
    26242653      // 1 bias increase by ones with slightly wrong djs
    26252654      // 2 bias by all
    2626       // 3 bias by all - tolerance
     2655      // 3 bias by all - tolerance 
    26272656#define TRYBIAS 3
    26282657
  • trunk/Test/CbcOrClpParam.cpp

    r735 r740  
    548548    case SPECIALOPTIONS:
    549549      model->setSpecialOptions(value);
     550    case THREADS:
     551      model->setNumberThreads(value);
    550552      break;
    551553    default:
     
    578580    value=model->specialOptions();
    579581    break;
     582  case THREADS:
     583    value = model->numberThreads();
    580584  default:
    581585    value=intValue_;
     
    22892293 elements may increase."
    22902294     );
     2295  parameters[numberParameters++]=
     2296    CbcOrClpParam("thread!s","Number of threads to try and use",
     2297                  -2,64,THREADS);
    22912298#endif
    22922299#ifdef COIN_USE_CBC
  • trunk/Test/CbcOrClpParam.hpp

    r731 r740  
    3434 
    3535  This coding scheme is in flux.
    36   NODESTRATEGY, BRANCHSTRATEGY, ADDCUTSSTRATEGY,
    37   CLEARCUTS, OSLSTUFF, CBCSTUFF are not used at present (03.10.24).
    3836*/
    3937
     
    5654    MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    5755    OUTPUTFORMAT,SLPVALUE,PRESOLVEOPTIONS,PRINTOPTIONS,SPECIALOPTIONS,
    58     SUBSTITUTION,DUALIZE,VERBOSE,
     56    SUBSTITUTION,DUALIZE,VERBOSE,THREADS,
    5957
    6058    STRONGBRANCHING=151,CUTDEPTH, MAXNODES,NUMBERBEFORE,NUMBERANALYZE,
  • trunk/include/ClpModel.hpp

    r722 r740  
    610610  inline void setWhatsChanged(int value)
    611611          { whatsChanged_ = value;} ;
     612  /// Number of threads (not really being used)
     613  inline int numberThreads() const
     614          { return numberThreads_;} ;
     615  inline void setNumberThreads(int value)
     616          { numberThreads_ = value;} ;
    612617  //@}
    613618  /**@name Message handling */
     
    883888  /// length of names (0 means no names)
    884889  int lengthNames_;
     890  /// Number of threads (not very operational)
     891  int numberThreads_;
    885892  /// Message handler
    886893  CoinMessageHandler * handler_;
  • trunk/include/ClpPackedMatrix.hpp

    r653 r740  
    1515    For details see CoinPackedMatrix */
    1616
     17class ClpPackedMatrix2;
    1718class ClpPackedMatrix : public ClpMatrixBase {
    1819 
     
    293294                    int numberRows, const int * whichRows,
    294295                    int numberColumns, const int * whichColumns) const ;
     296  /// make special row copy
     297  void specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy);
    295298   //@}
    296299   
     
    308311  /// Gaps flag - set true if column start and length don't say contiguous
    309312  bool hasGaps_;
     313  /// Special row copy
     314  ClpPackedMatrix2 * rowCopy_;
    310315   //@}
    311316};
    312 
     317#ifdef THREAD
     318#include <pthread.h>
     319typedef struct {
     320  double acceptablePivot;
     321  const ClpSimplex * model;
     322  double * spare;
     323  int * spareIndex;
     324  double * arrayTemp;
     325  int * indexTemp;
     326  int * numberInPtr;
     327  double * bestPossiblePtr;
     328  double * upperThetaPtr;
     329  int * posFreePtr;
     330  double * freePivotPtr;
     331  int * numberOutPtr;
     332  const unsigned short * count;
     333  const double * pi;
     334  const CoinBigIndex * rowStart;
     335  const double * element;
     336  const unsigned short * column;
     337  int offset;
     338  int numberInRowArray;
     339  int numberLook;
     340} dualColumn0Struct;
    313341#endif
     342class ClpPackedMatrix2 {
     343 
     344public:
     345  /**@name Useful methods */
     346  //@{
     347    /** Return <code>x * -1 * A in <code>z</code>.
     348        Note - x packed and z will be packed mode
     349        Squashes small elements and knows about ClpSimplex */
     350  void transposeTimes(const ClpSimplex * model,
     351                      const CoinPackedMatrix * rowCopy,
     352                      const CoinIndexedVector * x,
     353                      CoinIndexedVector * spareArray,
     354                      CoinIndexedVector * z) const;
     355  /// Returns true if copy has useful information
     356  inline bool usefulInfo() const
     357  { return rowStart_!=NULL;};
     358  //@}
     359
     360
     361  /**@name Constructors, destructor */
     362   //@{
     363   /** Default constructor. */
     364   ClpPackedMatrix2();
     365   /** Constructor from copy. */
     366   ClpPackedMatrix2(ClpSimplex * model,const CoinPackedMatrix * rowCopy);
     367   /** Destructor */
     368   virtual ~ClpPackedMatrix2();
     369   //@}
     370
     371   /**@name Copy method */
     372   //@{
     373   /** The copy constructor. */
     374   ClpPackedMatrix2(const ClpPackedMatrix2&);
     375   ClpPackedMatrix2& operator=(const ClpPackedMatrix2&);
     376   //@}
     377   
     378   
     379protected:
     380   /**@name Data members
     381      The data members are protected to allow access for derived classes. */
     382   //@{
     383  /// Number of blocks
     384  int numberBlocks_;
     385  /// Number of rows
     386  int numberRows_;
     387  /// Column offset for each block (plus one at end)
     388  int * offset_;
     389  /// Counts of elements in each part of row
     390  mutable unsigned short * count_;
     391  /// Row starts
     392  mutable CoinBigIndex * rowStart_;
     393  /// columns within block
     394  unsigned short * column_;
     395  /// work arrays
     396  double * work_;
     397#ifdef THREAD
     398  pthread_t * threadId_;
     399  dualColumn0Struct * info_;
     400#endif
     401   //@}
     402};
     403
     404#endif
  • trunk/include/ClpSimplex.hpp

    r737 r740  
    12641264  /// For dealing with all issues of cycling etc
    12651265  ClpSimplexProgress * progress_;
     1266public:
     1267  /// Spare int array for passing information [0]!=0 switches on
     1268  mutable int spareIntArray_[4];
     1269  /// Spare double array for passing information [0]!=0 switches on
     1270  mutable double spareDoubleArray_[4];
     1271protected:
    12661272  /// Allow OsiClp certain perks
    12671273  friend class OsiClpSolverInterface;
  • trunk/include/ClpSimplexDual.hpp

    r636 r740  
    196196                  double accpetablePivot,
    197197                  CoinBigIndex * dubiousWeights);
     198  /// Does first bit of dualColumn
     199  int dualColumn0(const CoinIndexedVector * rowArray,
     200                  const CoinIndexedVector * columnArray,
     201                  CoinIndexedVector * spareArray,
     202                  double acceptablePivot,
     203                  double & upperReturn, double &bestReturn);
    198204  /**
    199205      Row array has row part of pivot row
Note: See TracChangeset for help on using the changeset viewer.