Changeset 118 for trunk


Ignore:
Timestamp:
Feb 7, 2003 5:39:18 PM (17 years ago)
Author:
forrest
Message:

Adding Network matrix and PlusMinusOne?

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpFactorization.cpp

    r63 r118  
    77#include "ClpSimplex.hpp"
    88#include "ClpMatrixBase.hpp"
     9#include "ClpNetworkBasis.hpp"
    910
    1011
     
    1718//-------------------------------------------------------------------
    1819ClpFactorization::ClpFactorization () :
    19    CoinFactorization() {}
     20   CoinFactorization()
     21{
     22  networkBasis_ = NULL;
     23}
    2024
    2125//-------------------------------------------------------------------
     
    2327//-------------------------------------------------------------------
    2428ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
    25    CoinFactorization(rhs) {}
     29   CoinFactorization(rhs)
     30{
     31  if (rhs.networkBasis_)
     32    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
     33  else
     34    networkBasis_=NULL;
     35}
    2636
    2737ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
    28    CoinFactorization(rhs) {}
     38   CoinFactorization(rhs)
     39{
     40  networkBasis_=NULL;
     41}
    2942
    3043//-------------------------------------------------------------------
    3144// Destructor
    3245//-------------------------------------------------------------------
    33 ClpFactorization::~ClpFactorization () {}
     46ClpFactorization::~ClpFactorization ()
     47{
     48  delete networkBasis_;
     49}
    3450
    3551//----------------------------------------------------------------
     
    4157  if (this != &rhs) {
    4258    CoinFactorization::operator=(rhs);
     59    delete networkBasis_;
     60    if (rhs.networkBasis_)
     61      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
     62    else
     63      networkBasis_=NULL;
    4364  }
    4465  return *this;
     
    157178  return status_;
    158179}
     180/* Replaces one Column to basis,
     181   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
     182   If checkBeforeModifying is true will do all accuracy checks
     183   before modifying factorization.  Whether to set this depends on
     184   speed considerations.  You could just do this on first iteration
     185   after factorization and thereafter re-factorize
     186   partial update already in U */
     187int
     188ClpFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
     189                      int pivotRow,
     190                      double pivotCheck ,
     191                      bool checkBeforeModifying)
     192{
     193  if (!networkBasis_) {
     194    return CoinFactorization::replaceColumn(regionSparse,
     195                                            pivotRow,
     196                                            pivotCheck,
     197                                            checkBeforeModifying);
     198  } else {
     199    return networkBasis_->replaceColumn(regionSparse,
     200                                        pivotRow);
     201  }
     202}
     203
     204/* Updates one column (FTRAN) from region2
     205   number returned is negative if no room
     206   region1 starts as zero and is zero at end */
     207int
     208ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
     209                                 CoinIndexedVector * regionSparse2,
     210                                 bool FTUpdate)
     211{
     212  if (!networkBasis_) {
     213    return CoinFactorization::updateColumn(regionSparse,
     214                                           regionSparse2,
     215                                           FTUpdate);
     216  } else {
     217    return networkBasis_->updateColumn(regionSparse,
     218                                       regionSparse2);
     219  }
     220}
     221/* Updates one column (FTRAN) to/from array
     222   number returned is negative if no room
     223   ** For large problems you should ALWAYS know where the nonzeros
     224   are, so please try and migrate to previous method after you
     225   have got code working using this simple method - thank you!
     226   (the only exception is if you know input is dense e.g. rhs)
     227   region starts as zero and is zero at end */
     228int
     229ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
     230                        double array[] ) const
     231{
     232  if (!networkBasis_) {
     233    return CoinFactorization::updateColumn(regionSparse,
     234                                           array);
     235  } else {
     236    return networkBasis_->updateColumn(regionSparse,
     237                                                    array);
     238  }
     239}
     240/* Updates one column transpose (BTRAN)
     241   For large problems you should ALWAYS know where the nonzeros
     242   are, so please try and migrate to previous method after you
     243   have got code working using this simple method - thank you!
     244   (the only exception is if you know input is dense e.g. dense objective)
     245   returns number of nonzeros */
     246int
     247ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
     248                                          double array[] ) const
     249{
     250  if (!networkBasis_) {
     251    return CoinFactorization::updateColumnTranspose(regionSparse,
     252                                                    array);
     253  } else {
     254    return networkBasis_->updateColumnTranspose(regionSparse,
     255                                                    array);
     256  }
     257}
     258/* Updates one column (BTRAN) from region2
     259   region1 starts as zero and is zero at end */
     260int
     261ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
     262                          CoinIndexedVector * regionSparse2) const
     263{
     264  if (!networkBasis_) {
     265    return CoinFactorization::updateColumnTranspose(regionSparse,
     266                                                    regionSparse2);
     267  } else {
     268    return networkBasis_->updateColumnTranspose(regionSparse,
     269                                                    regionSparse2);
     270  }
     271}
  • trunk/ClpPackedMatrix.cpp

    r109 r118  
    220220  double zeroTolerance = model->factorization()->zeroTolerance();
    221221  int numberRows = model->numberRows();
    222   const ClpMatrixBase * rowCopy = model->rowCopy();
    223   // make sure row copy is of correct type
    224   // if not we would have to have another transposeTimes
    225   assert (!rowCopy||rowCopy->type()==1);
     222  ClpPackedMatrix* rowCopy =
     223    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
    226224  if (numberInRowArray>0.3*numberRows||!rowCopy) {
    227225    // do by column
     
    272270      }
    273271    }
    274   } else if (numberInRowArray>2||y->getNumElements()) {
     272    columnArray->setNumElements(numberNonZero);
     273    y->setNumElements(0);
     274  } else {
     275    // do by row
     276    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
     277  }
     278}
     279/* Return <code>x * A + y</code> in <code>z</code>.
     280        Squashes small elements and knows about ClpSimplex */
     281void
     282ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
     283                              const CoinIndexedVector * rowArray,
     284                              CoinIndexedVector * y,
     285                              CoinIndexedVector * columnArray) const
     286{
     287  columnArray->clear();
     288  double * pi = rowArray->denseVector();
     289  int numberNonZero=0;
     290  int * index = columnArray->getIndices();
     291  double * array = columnArray->denseVector();
     292  int numberInRowArray = rowArray->getNumElements();
     293  // maybe I need one in OsiSimplex
     294  double zeroTolerance = model->factorization()->zeroTolerance();
     295  const int * column = getIndices();
     296  const CoinBigIndex * rowStart = getVectorStarts();
     297  const double * element = getElements();
     298  const int * whichRow = rowArray->getIndices();
     299  if (numberInRowArray>2||y->getNumElements()) {
    275300    // do by rows
    276301    // ** Row copy is already scaled
     
    294319    }
    295320
    296     const int * column = rowCopy->getIndices();
    297     const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    298     const double * element = rowCopy->getElements();
    299     const int * whichRow = rowArray->getIndices();
    300321    for (i=0;i<numberInRowArray;i++) {
    301322      iRow = whichRow[i];
     
    330351    numberNonZero=0;
    331352
    332     const int * column = rowCopy->getIndices();
    333     const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    334     const double * element = rowCopy->getElements();
    335     const int * whichRow = rowArray->getIndices();
    336353    double value;
    337354    iRow = whichRow[0];
     
    371388    int iRow=rowArray->getIndices()[0];
    372389    numberNonZero=0;
    373     const int * column = rowCopy->getIndices();
    374     const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    375     const double * element = rowCopy->getElements();
    376390    double value = pi[iRow]*scalar;
    377391    CoinBigIndex j;
     
    712726  }
    713727}
    714 // Creates row copy and scales if necessary
    715 ClpMatrixBase *
    716 ClpPackedMatrix::scaledRowCopy(ClpSimplex * model) const
    717 {
    718   ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
    719   ClpPackedMatrix* rowCopy =
    720     dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
    721 
    722   // Make sure it is really a ClpPackedMatrix
    723   assert (rowCopy!=NULL);
    724 
    725   const double * rowScale = model->rowScale();
    726   const double * columnScale = model->columnScale();
    727 
    728   if (rowScale) {
    729     // scale row copy
    730     int numberRows = model->numberRows();
    731     int numberColumns = model->numberColumns();
    732     const int * column = rowCopy->getIndices();
    733     const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    734     const double * element = rowCopy->getElements();
    735     int iRow;
    736     // need to replace row by row
    737     double * newElement = new double[numberColumns];
    738     // scale row copy
    739     for (iRow=0;iRow<numberRows;iRow++) {
    740       int j;
    741       double scale = rowScale[iRow];
    742       const double * elementsInThisRow = element + rowStart[iRow];
    743       const int * columnsInThisRow = column + rowStart[iRow];
    744       int number = rowStart[iRow+1]-rowStart[iRow];
    745       assert (number<=numberColumns);
    746       for (j=0;j<number;j++) {
    747         int iColumn = columnsInThisRow[j];
    748         newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
    749       }
    750       rowCopy->replaceVector(iRow,number,newElement);
    751     }
    752     delete [] newElement;
    753   }
    754   return rowCopyBase;
    755 }
    756728/* Unpacks a column into an CoinIndexedvector
    757729      Note that model is NOT const.  Bounds and objective could
  • trunk/Makefile.Clp

    r109 r118  
    1818LIBSRC += ClpMessage.cpp
    1919LIBSRC += ClpModel.cpp
     20LIBSRC += ClpNetworkBasis.cpp
     21LIBSRC += ClpNetworkMatrix.cpp
    2022LIBSRC += ClpNonLinearCost.cpp
    2123LIBSRC += ClpPackedMatrix.cpp
     24LIBSRC += ClpPlusMinusOneMatrix.cpp
    2225LIBSRC += ClpPrimalColumnDantzig.cpp
    2326LIBSRC += ClpPrimalColumnPivot.cpp
  • trunk/Test/unitTest.cpp

    r115 r118  
    2424#include "ClpPrimalColumnDantzig.hpp"
    2525#include "ClpParameters.hpp"
     26#include "ClpNetworkMatrix.hpp"
     27#include "ClpPlusMinusOneMatrix.hpp"
    2628
    2729#include "Presolve.hpp"
     
    792794  }
    793795 
     796  // test network
     797  {   
     798    std::string fn = mpsDir+"input.130";
     799    int numberColumns;
     800    int numberRows;
     801   
     802    FILE * fp = fopen(fn.c_str(),"r");
     803
     804    if (!fp) {
     805      fprintf(stderr,"Unable to open file input.130 in mpsDir directory\n");
     806      exit(1);
     807    }
     808    int problem;
     809    char temp[100];
     810    // read and skip
     811    fscanf(fp,"%s",temp);
     812    assert (!strcmp(temp,"BEGIN"));
     813    fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows,
     814           &numberColumns);
     815    // scan down to SUPPLY
     816    while (fgets(temp,100,fp)) {
     817      if (!strncmp(temp,"SUPPLY",6))
     818        break;
     819    }
     820    if (strncmp(temp,"SUPPLY",6)) {
     821      fprintf(stderr,"Unable to find SUPPLY\n");
     822      exit(2);
     823    }
     824    // get space for rhs
     825    double * lower = new double[numberRows];
     826    double * upper = new double[numberRows];
     827    int i;
     828    for (i=0;i<numberRows;i++) {
     829      lower[i]=0.0;
     830      upper[i]=0.0;
     831    }
     832    // ***** Remember to convert to C notation
     833    while (fgets(temp,100,fp)) {
     834      int row;
     835      int value;
     836      if (!strncmp(temp,"ARCS",4))
     837        break;
     838      sscanf(temp,"%d %d",&row,&value);
     839      upper[row-1]=-value;
     840      lower[row-1]=-value;
     841    }
     842    if (strncmp(temp,"ARCS",4)) {
     843      fprintf(stderr,"Unable to find ARCS\n");
     844      exit(2);
     845    }
     846    // number of columns may be underestimate
     847    int * head = new int[2*numberColumns];
     848    int * tail = new int[2*numberColumns];
     849    int * ub = new int[2*numberColumns];
     850    int * cost = new int[2*numberColumns];
     851    // ***** Remember to convert to C notation
     852    numberColumns=0;
     853    while (fgets(temp,100,fp)) {
     854      int iHead;
     855      int iTail;
     856      int iUb;
     857      int iCost;
     858      if (!strncmp(temp,"DEMAND",6))
     859        break;
     860      sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
     861      iHead--;
     862      iTail--;
     863      head[numberColumns]=iHead;
     864      tail[numberColumns]=iTail;
     865      ub[numberColumns]=iUb;
     866      cost[numberColumns]=iCost;
     867      numberColumns++;
     868    }
     869    if (strncmp(temp,"DEMAND",6)) {
     870      fprintf(stderr,"Unable to find DEMAND\n");
     871      exit(2);
     872    }
     873    // ***** Remember to convert to C notation
     874    while (fgets(temp,100,fp)) {
     875      int row;
     876      int value;
     877      if (!strncmp(temp,"END",3))
     878        break;
     879      sscanf(temp,"%d %d",&row,&value);
     880      upper[row-1]=value;
     881      lower[row-1]=value;
     882    }
     883    if (strncmp(temp,"END",3)) {
     884      fprintf(stderr,"Unable to find END\n");
     885      exit(2);
     886    }
     887    printf("Problem %d has %d rows and %d columns\n",problem,
     888           numberRows,numberColumns);
     889    fclose(fp);
     890    ClpSimplex  model;
     891    // now build model
     892   
     893    double * objective =new double[numberColumns];
     894    double * lowerColumn = new double[numberColumns];
     895    double * upperColumn = new double[numberColumns];
     896   
     897    double * element = new double [2*numberColumns];
     898    int * start = new int[numberColumns+1];
     899    int * row = new int[2*numberColumns];
     900    start[numberColumns]=2*numberColumns;
     901    for (i=0;i<numberColumns;i++) {
     902      start[i]=2*i;
     903      element[2*i]=-1.0;
     904      element[2*i+1]=1.0;
     905      row[2*i]=head[i];
     906      row[2*i+1]=tail[i];
     907      lowerColumn[i]=0.0;
     908      upperColumn[i]=ub[i];
     909      objective[i]=cost[i];
     910    }
     911    // Create Packed Matrix
     912    CoinPackedMatrix matrix;
     913    int * lengths = NULL;
     914    matrix.assignMatrix(true,numberRows,numberColumns,
     915                        2*numberColumns,element,row,start,lengths);
     916    // load model
     917   
     918    model.loadProblem(matrix,
     919                      lowerColumn,upperColumn,objective,
     920                      lower,upper);
     921   
     922    model.factorization()->maximumPivots(200+model.numberRows()/100);
     923    model.createStatus();
     924    double time1 = cpuTime();
     925    model.dual();
     926    std::cout<<"Network problem, ClpPackedMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
     927    ClpPlusMinusOneMatrix plusMinus(matrix);
     928    assert (plusMinus.getIndices()); // would be zero if not +- one
     929    model.loadProblem(plusMinus,
     930                      lowerColumn,upperColumn,objective,
     931                      lower,upper);
     932   
     933    model.factorization()->maximumPivots(200+model.numberRows()/100);
     934    model.createStatus();
     935    time1 = cpuTime();
     936    model.dual();
     937    std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
     938    ClpNetworkMatrix network(numberColumns,head,tail);
     939    model.loadProblem(network,
     940                      lowerColumn,upperColumn,objective,
     941                      lower,upper);
     942   
     943    model.factorization()->maximumPivots(200+model.numberRows()/100);
     944    model.createStatus();
     945    time1 = cpuTime();
     946    model.dual();
     947    std::cout<<"Network problem, ClpNetworkMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
     948    delete [] lower;
     949    delete [] upper;
     950    delete [] head;
     951    delete [] tail;
     952    delete [] ub;
     953    delete [] cost;
     954    delete [] objective;
     955    delete [] lowerColumn;
     956    delete [] upperColumn;
     957  }
     958 
    794959}
  • trunk/include/ClpFactorization.hpp

    r63 r118  
    1010
    1111/** This just implements CoinFactorization when an ClpMatrixBase object
    12     is passed.  It has no data.
     12    is passed.  If a network then has a dummy CoinFactorization and
     13    a genuine ClpNetworkBasis object
    1314*/
    1415class ClpMatrixBase;
    1516class ClpSimplex;
     17class ClpNetworkBasis;
    1618
    1719class ClpFactorization : public CoinFactorization {
     
    5557   //@}
    5658   
     59  /*  **** below here is so can use networkish basis */
     60  /**@name rank one updates which do exist */
     61  //@{
     62
     63  /** Replaces one Column to basis,
     64   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
     65      If checkBeforeModifying is true will do all accuracy checks
     66      before modifying factorization.  Whether to set this depends on
     67      speed considerations.  You could just do this on first iteration
     68      after factorization and thereafter re-factorize
     69   partial update already in U */
     70  int replaceColumn ( CoinIndexedVector * regionSparse,
     71                      int pivotRow,
     72                      double pivotCheck ,
     73                      bool checkBeforeModifying=false);
     74  //@}
     75
     76  /**@name various uses of factorization (return code number elements)
     77   which user may want to know about */
     78  //@{
     79  /** Updates one column (FTRAN) from region2
     80      number returned is negative if no room
     81      region1 starts as zero and is zero at end */
     82  int updateColumn ( CoinIndexedVector * regionSparse,
     83                        CoinIndexedVector * regionSparse2,
     84                        bool FTUpdate = false ) ;
     85  /** Updates one column (FTRAN) to/from array
     86      number returned is negative if no room
     87      ** For large problems you should ALWAYS know where the nonzeros
     88      are, so please try and migrate to previous method after you
     89      have got code working using this simple method - thank you!
     90      (the only exception is if you know input is dense e.g. rhs)
     91      region starts as zero and is zero at end */
     92  int updateColumn ( CoinIndexedVector * regionSparse,
     93                        double array[] ) const;
     94  /** Updates one column transpose (BTRAN)
     95      ** For large problems you should ALWAYS know where the nonzeros
     96      are, so please try and migrate to previous method after you
     97      have got code working using this simple method - thank you!
     98      (the only exception is if you know input is dense e.g. dense objective)
     99      returns number of nonzeros */
     100  int updateColumnTranspose ( CoinIndexedVector * regionSparse,
     101                                 double array[] ) const;
     102  /** Updates one column (BTRAN) from region2
     103      region1 starts as zero and is zero at end */
     104  int updateColumnTranspose ( CoinIndexedVector * regionSparse,
     105                              CoinIndexedVector * regionSparse2) const;
     106  //@}
    57107   
     108////////////////// data //////////////////
     109private:
     110
     111  /**@name data */
     112  //@{
     113  /// Pointer to network basis
     114  ClpNetworkBasis * networkBasis_;
     115  //@}
    58116};
    59117
  • trunk/include/ClpMatrixBase.hpp

    r78 r118  
    7676  virtual int scale(ClpSimplex * model) const
    7777  { return 1;};
    78   /// Creates row copy and scales if necessary
    79   virtual ClpMatrixBase * scaledRowCopy(ClpSimplex * model) const
    80   { return reverseOrderedCopy();};
    8178
    8279  /** Checks if all elements are in valid range.  Can just
  • trunk/include/ClpPackedMatrix.hpp

    r78 r118  
    7878      returns non-zero if no scaling done */
    7979  virtual int scale(ClpSimplex * model) const ;
    80   /// Creates row copy and scales if necessary
    81   virtual ClpMatrixBase * scaledRowCopy(ClpSimplex * model) const;
    8280  /** Checks if all elements are in valid range.  Can just
    8381      return true if you are not paranoid.  For Clp I will
     
    127125        Squashes small elements and knows about ClpSimplex */
    128126  virtual void transposeTimes(const ClpSimplex * model, double scalar,
     127                              const CoinIndexedVector * x,
     128                              CoinIndexedVector * y,
     129                              CoinIndexedVector * z) const;
     130    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
     131        Can use y as temporary array (will be empty at end)
     132        Squashes small elements and knows about ClpSimplex.
     133    This version uses row copy*/
     134  virtual void transposeTimesByRow(const ClpSimplex * model, double scalar,
    129135                              const CoinIndexedVector * x,
    130136                              CoinIndexedVector * y,
Note: See TracChangeset for help on using the changeset viewer.