Ignore:
Timestamp:
Oct 15, 2003 5:34:57 PM (16 years ago)
Author:
forrest
Message:

This should break everything

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpNetworkMatrix.cpp

    r124 r225  
    1515#include "ClpPlusMinusOneMatrix.hpp"
    1616#include "ClpMessage.hpp"
     17#include <iostream>
    1718
    1819//#############################################################################
     
    369370  ClpPlusMinusOneMatrix* rowCopy =
    370371    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    371   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     372  bool packed = rowArray->packedMode();
     373  double factor = 0.3;
     374  // We may not want to do by row if there may be cache problems
     375  int numberColumns = model->numberColumns();
     376  // It would be nice to find L2 cache size - for moment 512K
     377  // Be slightly optimistic
     378  if (numberColumns*sizeof(double)>1000000) {
     379    if (numberRows*10<numberColumns)
     380      factor=0.1;
     381    else if (numberRows*4<numberColumns)
     382      factor=0.15;
     383    else if (numberRows*2<numberColumns)
     384      factor=0.2;
     385    //if (model->numberIterations()%50==0)
     386    //printf("%d nonzero\n",numberInRowArray);
     387  }
     388  if (numberInRowArray>factor*numberRows||!rowCopy) {
    372389    // do by column
    373390    int iColumn;
    374     double * markVector = y->denseVector(); // probably empty
     391    assert (!y->getNumElements());
    375392    CoinBigIndex j=0;
    376     if (trueNetwork_) {
    377       for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    378         double value = markVector[iColumn];
    379         markVector[iColumn]=0.0;
    380         int iRowM = indices_[j];
    381         int iRowP = indices_[j+1];
    382         value -= scalar*pi[iRowM];
    383         value += scalar*pi[iRowP];
    384         if (fabs(value)>zeroTolerance) {
    385           index[numberNonZero++]=iColumn;
    386           array[iColumn]=value;
    387         }
     393    if (packed) {
     394      // need to expand pi into y
     395      assert(y->capacity()>=numberRows);
     396      double * piOld = pi;
     397      pi = y->denseVector();
     398      const int * whichRow = rowArray->getIndices();
     399      int i;
     400      // modify pi so can collapse to one loop
     401      for (i=0;i<numberInRowArray;i++) {
     402        int iRow = whichRow[i];
     403        pi[iRow]=scalar*piOld[i];
     404      }
     405      if (trueNetwork_) {
     406        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     407          double value = 0.0;
     408          int iRowM = indices_[j];
     409          int iRowP = indices_[j+1];
     410          value -= pi[iRowM];
     411          value += pi[iRowP];
     412          if (fabs(value)>zeroTolerance) {
     413            array[numberNonZero]=value;
     414            index[numberNonZero++]=iColumn;
     415          }
     416        }
     417      } else {
     418        // skip negative rows
     419        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     420          double value = 0.0;
     421          int iRowM = indices_[j];
     422          int iRowP = indices_[j+1];
     423          if (iRowM>=0)
     424            value -= pi[iRowM];
     425          if (iRowP>=0)
     426            value += pi[iRowP];
     427          if (fabs(value)>zeroTolerance) {
     428            array[numberNonZero]=value;
     429            index[numberNonZero++]=iColumn;
     430          }
     431        }
     432      }
     433      for (i=0;i<numberInRowArray;i++) {
     434        int iRow = whichRow[i];
     435        pi[iRow]=0.0;
    388436      }
    389437    } else {
    390       // skip negative rows
    391       for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    392         double value = markVector[iColumn];
    393         markVector[iColumn]=0.0;
    394         int iRowM = indices_[j];
    395         int iRowP = indices_[j+1];
    396         if (iRowM>=0)
     438      if (trueNetwork_) {
     439        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     440          double value = 0.0;
     441          int iRowM = indices_[j];
     442          int iRowP = indices_[j+1];
    397443          value -= scalar*pi[iRowM];
    398         if (iRowP>=0)
    399444          value += scalar*pi[iRowP];
    400         if (fabs(value)>zeroTolerance) {
    401           index[numberNonZero++]=iColumn;
    402           array[iColumn]=value;
     445          if (fabs(value)>zeroTolerance) {
     446            index[numberNonZero++]=iColumn;
     447            array[iColumn]=value;
     448          }
     449        }
     450      } else {
     451        // skip negative rows
     452        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     453          double value = 0.0;
     454          int iRowM = indices_[j];
     455          int iRowP = indices_[j+1];
     456          if (iRowM>=0)
     457            value -= scalar*pi[iRowM];
     458          if (iRowP>=0)
     459            value += scalar*pi[iRowP];
     460          if (fabs(value)>zeroTolerance) {
     461            index[numberNonZero++]=iColumn;
     462            array[iColumn]=value;
     463          }
    403464        }
    404465      }
    405466    }
    406467    columnArray->setNumElements(numberNonZero);
    407     y->setNumElements(0);
    408468  } else {
    409469    // do by row
     
    430490  int numberToDo = y->getNumElements();
    431491  const int * which = y->getIndices();
    432   if (trueNetwork_) {
    433     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    434       int iColumn = which[jColumn];
    435       double value = 0.0;
    436       CoinBigIndex j=iColumn<<1;
    437       int iRowM = indices_[j];
    438       int iRowP = indices_[j+1];
    439       value -= pi[iRowM];
    440       value += pi[iRowP];
    441       if (fabs(value)>zeroTolerance) {
    442         index[numberNonZero++]=iColumn;
    443         array[iColumn]=value;
    444       }
     492  bool packed = rowArray->packedMode();
     493  if (packed) {
     494    // Must line up with y
     495    // need to expand pi into y
     496    int numberInRowArray = rowArray->getNumElements();
     497    assert(y->capacity()>=model->numberRows());
     498    double * piOld = pi;
     499    pi = y->denseVector();
     500    const int * whichRow = rowArray->getIndices();
     501    int i;
     502    for (i=0;i<numberInRowArray;i++) {
     503      int iRow = whichRow[i];
     504      pi[iRow]=piOld[i];
     505    }
     506    if (trueNetwork_) {
     507      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     508        int iColumn = which[jColumn];
     509        double value = 0.0;
     510        CoinBigIndex j=iColumn<<1;
     511        int iRowM = indices_[j];
     512        int iRowP = indices_[j+1];
     513        value -= pi[iRowM];
     514        value += pi[iRowP];
     515        array[jColumn]=value;
     516      }
     517    } else {
     518      // skip negative rows
     519      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     520        int iColumn = which[jColumn];
     521        double value = 0.0;
     522        CoinBigIndex j=iColumn<<1;
     523        int iRowM = indices_[j];
     524        int iRowP = indices_[j+1];
     525        if (iRowM>=0)
     526          value -= pi[iRowM];
     527        if (iRowP>=0)
     528          value += pi[iRowP];
     529        array[jColumn]=value;
     530      }
     531    }
     532    for (i=0;i<numberInRowArray;i++) {
     533      int iRow = whichRow[i];
     534      pi[iRow]=0.0;
    445535    }
    446536  } else {
    447     // skip negative rows
    448     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    449       int iColumn = which[jColumn];
    450       double value = 0.0;
    451       CoinBigIndex j=iColumn<<1;
    452       int iRowM = indices_[j];
    453       int iRowP = indices_[j+1];
    454       if (iRowM>=0)
     537    if (trueNetwork_) {
     538      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     539        int iColumn = which[jColumn];
     540        double value = 0.0;
     541        CoinBigIndex j=iColumn<<1;
     542        int iRowM = indices_[j];
     543        int iRowP = indices_[j+1];
    455544        value -= pi[iRowM];
    456       if (iRowP>=0)
    457545        value += pi[iRowP];
    458       if (fabs(value)>zeroTolerance) {
    459         index[numberNonZero++]=iColumn;
    460         array[iColumn]=value;
     546        if (fabs(value)>zeroTolerance) {
     547          index[numberNonZero++]=iColumn;
     548          array[iColumn]=value;
     549        }
     550      }
     551    } else {
     552      // skip negative rows
     553      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     554        int iColumn = which[jColumn];
     555        double value = 0.0;
     556        CoinBigIndex j=iColumn<<1;
     557        int iRowM = indices_[j];
     558        int iRowP = indices_[j+1];
     559        if (iRowM>=0)
     560          value -= pi[iRowM];
     561        if (iRowP>=0)
     562          value += pi[iRowP];
     563        if (fabs(value)>zeroTolerance) {
     564          index[numberNonZero++]=iColumn;
     565          array[iColumn]=value;
     566        }
    461567      }
    462568    }
     
    542648  return numberElements;
    543649}
     650/* If element NULL returns number of elements in column part of basis,
     651   If not NULL fills in as well */
     652CoinBigIndex
     653ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
     654                                 const int * whichColumn,
     655                                 int numberBasic,
     656                                 int numberColumnBasic,
     657                                 int * indexRowU, int * indexColumnU,
     658                                 double * elementU) const
     659{
     660  int i;
     661  CoinBigIndex numberElements=0;
     662  if (elementU!=NULL) {
     663    if (trueNetwork_) {
     664      for (i=0;i<numberColumnBasic;i++) {
     665        int iColumn = whichColumn[i];
     666        CoinBigIndex j=iColumn<<1;
     667        int iRowM = indices_[j];
     668        int iRowP = indices_[j+1];
     669        indexRowU[numberElements]=iRowM;
     670        indexColumnU[numberElements]=numberBasic;
     671        elementU[numberElements]=-1.0;
     672        indexRowU[numberElements+1]=iRowP;
     673        indexColumnU[numberElements+1]=numberBasic;
     674        elementU[numberElements+1]=1.0;
     675        numberElements+=2;
     676        numberBasic++;
     677      }
     678    } else {
     679      for (i=0;i<numberColumnBasic;i++) {
     680        int iColumn = whichColumn[i];
     681        CoinBigIndex j=iColumn<<1;
     682        int iRowM = indices_[j];
     683        int iRowP = indices_[j+1];
     684        if (iRowM>=0) {
     685          indexRowU[numberElements]=iRowM;
     686          indexColumnU[numberElements]=numberBasic;
     687          elementU[numberElements++]=-1.0;
     688        }
     689        if (iRowP>=0) {
     690          indexRowU[numberElements]=iRowP;
     691          indexColumnU[numberElements]=numberBasic;
     692          elementU[numberElements++]=1.0;
     693        }
     694        numberBasic++;
     695      }
     696    }
     697  } else {
     698    if (trueNetwork_) {
     699      numberElements = 2*numberColumnBasic;
     700    } else {
     701      for (i=0;i<numberColumnBasic;i++) {
     702        int iColumn = whichColumn[i];
     703        CoinBigIndex j=iColumn<<1;
     704        int iRowM = indices_[j];
     705        int iRowP = indices_[j+1];
     706        if (iRowM>=0)
     707          numberElements ++;
     708        if (iRowP>=0)
     709          numberElements ++;
     710      }
     711    }
     712  }
     713  return numberElements;
     714}
    544715/* Unpacks a column into an CoinIndexedvector
    545       Note that model is NOT const.  Bounds and objective could
    546       be modified if doing column generation */
     716 */
    547717void
    548718ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
     
    557727    rowArray->add(iRowP,1.0);
    558728}
     729/* Unpacks a column into an CoinIndexedvector
     730** in packed foramt
     731Note that model is NOT const.  Bounds and objective could
     732be modified if doing column generation (just for this variable) */
     733void
     734ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
     735                            CoinIndexedVector * rowArray,
     736                            int iColumn) const
     737{
     738  int * index = rowArray->getIndices();
     739  double * array = rowArray->denseVector();
     740  int number = 0;
     741  CoinBigIndex j=iColumn<<1;
     742  int iRowM = indices_[j];
     743  int iRowP = indices_[j+1];
     744  if (iRowM>=0) {
     745    array[number]=-1.0;
     746    index[number++]=iRowM;
     747  }
     748  if (iRowP>=0) {
     749    array[number]=1.0;
     750    index[number++]=iRowP;
     751  }
     752  rowArray->setNumElements(number);
     753  rowArray->setPackedMode(true);
     754}
    559755/* Adds multiple of a column into an CoinIndexedvector
    560756      You can use quickAdd to add to vector */
     
    631827void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel)
    632828{
     829  std::cerr<<"deleteCols not implemented in ClpNetworkMatrix"<<std::endl;
    633830  abort();
    634831}
     
    636833void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel)
    637834{
     835  std::cerr<<"deleteRows not implemented in ClpNetworkMatrix"<<std::endl;
    638836  abort();
    639837}
     838/* Given positive integer weights for each row fills in sum of weights
     839   for each column (and slack).
     840   Returns weights vector
     841*/
     842CoinBigIndex *
     843ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
     844{
     845  int numberRows = model->numberRows();
     846  int numberColumns =model->numberColumns();
     847  int number = numberRows+numberColumns;
     848  CoinBigIndex * weights = new CoinBigIndex[number];
     849  int i;
     850  for (i=0;i<numberColumns;i++) {
     851    CoinBigIndex j=i<<1;
     852    CoinBigIndex count=0;
     853    int iRowM = indices_[j];
     854    int iRowP = indices_[j+1];
     855    if (iRowM>=0) {
     856      count += inputWeights[iRowM];
     857    }
     858    if (iRowP>=0) {
     859      count += inputWeights[iRowP];
     860    }
     861    weights[i]=count;
     862  }
     863  for (i=0;i<numberRows;i++) {
     864    weights[i+numberColumns]=inputWeights[i];
     865  }
     866  return weights;
     867}
     868/* Returns largest and smallest elements of both signs.
     869   Largest refers to largest absolute value.
     870*/
     871void
     872ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
     873                       double & smallestPositive, double & largestPositive)
     874{
     875  smallestNegative=-1.0;
     876  largestNegative=-1.0;
     877  smallestPositive=1.0;
     878  largestPositive=1.0;
     879}
Note: See TracChangeset for help on using the changeset viewer.