Changeset 119


Ignore:
Timestamp:
Feb 9, 2003 6:12:40 PM (17 years ago)
Author:
forrest
Message:

Updating Clp stuff for networks

Files:
19 edited

Legend:

Unmodified
Added
Removed
  • branches/unlabeled-1.1.2/ClpNetworkBasis.cpp

    r118 r119  
    77#include "ClpSimplex.hpp"
    88#include "ClpMatrixBase.hpp"
     9#include "CoinIndexedVector.hpp"
    910
    1011
     
    1819ClpNetworkBasis::ClpNetworkBasis ()
    1920{
     21  slackValue_=-1.0;
     22  numberRows_=0;
     23  numberColumns_=0;
     24  root_ = -1;
     25  leaf_ = -1;
     26  parent_ = NULL;
     27  descendant_ = NULL;
     28  pivot_ = NULL;
     29  rightSibling_ = NULL;
     30  leftSibling_ = NULL;
     31  sign_ = NULL;
     32  stack_ = NULL;
     33  toLeaf_ = NULL;
     34  toRoot_ = NULL;
     35  mark_ = NULL;
     36  model_=NULL;
     37}
     38// Constructor from CoinFactorization
     39ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex * model,
     40                                 int numberRows, const double * pivotRegion,
     41                                 const int * permuteBack,
     42                                 const int * startColumn,
     43                                 const int * numberInColumn,
     44                                 const int * indexRow, const double * element)
     45{
     46  slackValue_=-1.0;
     47  numberRows_=numberRows;
     48  numberColumns_=numberRows;
     49  parent_ = new int [ numberRows_+1];
     50  descendant_ = new int [ numberRows_+1];
     51  pivot_ = new int [ numberRows_+1];
     52  rightSibling_ = new int [ numberRows_+1];
     53  leftSibling_ = new int [ numberRows_+1];
     54  sign_ = new double [ numberRows_+1];
     55  stack_ = new int [ numberRows_+1];
     56  toLeaf_ = new int [numberRows_+1];
     57  toRoot_ = new int [numberRows_+1];
     58  mark_ = new char[numberRows_+1];
     59  int i;
     60  for (i=0;i<numberRows_+1;i++) {
     61    parent_[i]=-1;
     62    descendant_[i]=-1;
     63    pivot_[i]=-1;
     64    rightSibling_[i]=-1;
     65    leftSibling_[i]=-1;
     66    sign_[i]=-1.0;
     67    stack_[i]=-1;
     68    mark_[i]=0;
     69  }
     70  // pivotColumnBack gives order of pivoting into basis
     71  // so pivotColumnback[0] is first slack in basis and
     72  // it pivots on row permuteBack[0]
     73  // a known root is given by permuteBack[numberRows_-1]
     74  root_ = numberRows_;
     75  int lastPivot=numberRows_;
     76  for (i=0;i<numberRows_;i++) {
     77    int iPivot = permuteBack[i];
     78    toRoot_[iPivot] = lastPivot;
     79    toLeaf_[lastPivot]=iPivot;
     80    lastPivot=iPivot;
     81    double sign;
     82    if (pivotRegion[i]>0.0)
     83      sign = 1.0;
     84    else
     85      sign =-1.0;
     86    int other;
     87    if (numberInColumn[i]>0) {
     88      int iRow = indexRow[startColumn[i]];
     89      other = permuteBack[iRow];
     90      assert (parent_[other]!=-1);
     91    } else {
     92      other = numberRows_;
     93    }
     94    sign_[iPivot] = sign;
     95    int iParent = other;
     96    parent_[iPivot] = other;
     97    if (descendant_[iParent]>=0) {
     98      // we have a sibling
     99      int iRight = descendant_[iParent];
     100      rightSibling_[iPivot]=iRight;
     101      leftSibling_[iRight]=iPivot;
     102    } else {
     103      rightSibling_[iPivot]=-1;
     104    }       
     105    descendant_[iParent] = iPivot;
     106    leftSibling_[iPivot]=-1;
     107  }
     108  toLeaf_[lastPivot]=numberRows_;
     109  leaf_ = lastPivot;
     110  model_=model;
    20111}
    21112
     
    25116ClpNetworkBasis::ClpNetworkBasis (const ClpNetworkBasis & rhs)
    26117{
     118  slackValue_=rhs.slackValue_;
     119  numberRows_=rhs.numberRows_;
     120  numberColumns_=rhs.numberColumns_;
     121  root_ = rhs.root_;
     122  leaf_ = rhs.leaf_;
     123  if (rhs.parent_) {
     124    parent_ = new int [numberRows_+1];
     125    memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
     126  } else {
     127    parent_ = NULL;
     128  }
     129  if (rhs.descendant_) {
     130    descendant_ = new int [numberRows_+1];
     131    memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
     132  } else {
     133    descendant_ = NULL;
     134  }
     135  if (rhs.pivot_) {
     136    pivot_ = new int [numberRows_+1];
     137    memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
     138  } else {
     139    pivot_ = NULL;
     140  }
     141  if (rhs.rightSibling_) {
     142    rightSibling_ = new int [numberRows_+1];
     143    memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
     144  } else {
     145    rightSibling_ = NULL;
     146  }
     147  if (rhs.leftSibling_) {
     148    leftSibling_ = new int [numberRows_+1];
     149    memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
     150  } else {
     151    leftSibling_ = NULL;
     152  }
     153  if (rhs.sign_) {
     154    sign_ = new double [numberRows_+1];
     155    memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
     156  } else {
     157    sign_ = NULL;
     158  }
     159  if (rhs.stack_) {
     160    stack_ = new int [numberRows_+1];
     161    memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
     162  } else {
     163    stack_ = NULL;
     164  }
     165  if (rhs.toLeaf_) {
     166    toLeaf_ = new int [numberRows_+1];
     167    memcpy(toLeaf_,rhs.toLeaf_,(numberRows_+1)*sizeof(int));
     168  } else {
     169    toLeaf_ = NULL;
     170  }
     171  if (rhs.toRoot_) {
     172    toRoot_ = new int [numberRows_+1];
     173    memcpy(toRoot_,rhs.toRoot_,(numberRows_+1)*sizeof(int));
     174  } else {
     175    toRoot_ = NULL;
     176  }
     177  if (rhs.mark_) {
     178    mark_ = new char [numberRows_+1];
     179    memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
     180  } else {
     181    mark_ = NULL;
     182  }
     183  model_=rhs.model_;
    27184}
    28185
     
    32189ClpNetworkBasis::~ClpNetworkBasis ()
    33190{
     191  delete [] parent_;
     192  delete [] descendant_;
     193  delete [] pivot_;
     194  delete [] rightSibling_;
     195  delete [] leftSibling_;
     196  delete [] sign_;
     197  delete [] stack_;
     198  delete [] toLeaf_;
     199  delete [] toRoot_;
     200  delete [] mark_;
    34201}
    35202
     
    41208{
    42209  if (this != &rhs) {
     210    delete [] parent_;
     211    delete [] descendant_;
     212    delete [] pivot_;
     213    delete [] rightSibling_;
     214    delete [] leftSibling_;
     215    delete [] sign_;
     216    delete [] stack_;
     217    delete [] toLeaf_;
     218    delete [] toRoot_;
     219    delete [] mark_;
     220    slackValue_=rhs.slackValue_;
     221    numberRows_=rhs.numberRows_;
     222    numberColumns_=rhs.numberColumns_;
     223    root_ = rhs.root_;
     224    leaf_ = rhs.leaf_;
     225    if (rhs.parent_) {
     226      parent_ = new int [numberRows_+1];
     227      memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
     228    } else {
     229      parent_ = NULL;
     230    }
     231    if (rhs.descendant_) {
     232      descendant_ = new int [numberRows_+1];
     233      memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
     234    } else {
     235      descendant_ = NULL;
     236    }
     237    if (rhs.pivot_) {
     238      pivot_ = new int [numberRows_+1];
     239      memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
     240    } else {
     241      pivot_ = NULL;
     242    }
     243    if (rhs.rightSibling_) {
     244      rightSibling_ = new int [numberRows_+1];
     245      memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
     246    } else {
     247      rightSibling_ = NULL;
     248    }
     249    if (rhs.leftSibling_) {
     250      leftSibling_ = new int [numberRows_+1];
     251      memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
     252    } else {
     253      leftSibling_ = NULL;
     254    }
     255    if (rhs.sign_) {
     256      sign_ = new double [numberRows_+1];
     257      memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
     258    } else {
     259      sign_ = NULL;
     260    }
     261    if (rhs.stack_) {
     262      stack_ = new int [numberRows_+1];
     263      memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
     264    } else {
     265      stack_ = NULL;
     266    }
     267    if (rhs.toLeaf_) {
     268      toLeaf_ = new int [numberRows_+1];
     269      memcpy(toLeaf_,rhs.toLeaf_,(numberRows_+1)*sizeof(int));
     270    } else {
     271      toLeaf_ = NULL;
     272    }
     273    if (rhs.toRoot_) {
     274      toRoot_ = new int [numberRows_+1];
     275      memcpy(toRoot_,rhs.toRoot_,(numberRows_+1)*sizeof(int));
     276    } else {
     277      toRoot_ = NULL;
     278    }
     279    if (rhs.mark_) {
     280      mark_ = new char [numberRows_+1];
     281      memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
     282    } else {
     283    mark_ = NULL;
     284    }
     285    model_=rhs.model_;
    43286  }
    44287  return *this;
    45288}
    46289/* Replaces one Column to basis,
    47    returns 0=OK, 1=Probably OK, 2=singular
     290   returns 0=OK
    48291*/
    49292int
     
    51294                                 int pivotRow)
    52295{
    53   abort();
    54   return 1;
    55 }
    56 
    57 /* Updates one column (FTRAN) from region2
    58    number returned is negative if no room
    59    region1 starts as zero and is zero at end */
     296  // regionSparse is empty
     297  assert (!regionSparse->getNumElements());
     298  model_->unpack(regionSparse, model_->sequenceIn());
     299  // arc given by pivotRow is leaving basis
     300  int iParent = parent_[pivotRow];
     301  // arc coming in has these two nodes
     302  int * indices = regionSparse->getIndices();
     303  int iRow0 = indices[0];
     304  int iRow1;
     305  if (regionSparse->getNumElements()==2)
     306    iRow1 = indices[1];
     307  else
     308    iRow1 = numberRows_;
     309  double sign = -regionSparse->denseVector()[iRow0];
     310  printf("In %d (%g) %d pivoting on %d\n",
     311         iRow1, sign, iRow0,pivotRow);
     312  regionSparse->clear();
     313  // take out of tree
     314  int iLeft = leftSibling_[pivotRow];
     315  int iRight = rightSibling_[pivotRow];
     316  if (iLeft>=0) {
     317    rightSibling_[iLeft] = iRight;
     318    if (iRight>=0)
     319      leftSibling_[iRight]=iLeft;
     320  } else if (iRight>=0) {
     321    leftSibling_[iRight]=-1;
     322    descendant_[iParent]=iRight;
     323  } else {
     324    descendant_[iParent]=-1;;
     325  }
     326  // move to other end of chain
     327  int descendant = descendant_[pivotRow];
     328  if (descendant>=0) {
     329    // make this descendant of that
     330    if (descendant_[descendant]>=0) {
     331      // we have a sibling
     332      int iRight = descendant_[descendant];
     333      rightSibling_[pivotRow]=iRight;
     334      leftSibling_[iRight]=pivotRow;
     335    } else {
     336      rightSibling_[pivotRow]=-1;
     337    }       
     338    descendant_[descendant] = pivotRow;
     339    leftSibling_[pivotRow]=-1;
     340  }
     341  // now insert new one
     342  descendant = descendant_[iRow1];
     343  parent_[iRow0] = iRow1;
     344  sign_[iRow1]= sign;
     345  if (descendant>=0) {
     346    // we have a sibling
     347    int iRight = descendant;
     348    rightSibling_[iRow1]=iRight;
     349    leftSibling_[iRight]=iRow1;
     350  } else {
     351    rightSibling_[iRow1]=-1;
     352  }         
     353  descendant_[descendant] = iRow1;
     354  leftSibling_[iRow1]=-1;
     355  return 0;
     356}
     357
     358/* Updates one column (FTRAN) from region2 */
    60359int
    61 ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse,
    62                                 CoinIndexedVector * regionSparse2)
    63 {
    64   abort();
    65   return 1;
     360ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse2)
     361{
     362  int iPivot=leaf_;
     363  int numberNonZero=0;
     364  double * array = regionSparse2->denseVector();
     365  while (iPivot!=numberRows_) {
     366    double pivotValue = array[iPivot];
     367    if (pivotValue) {
     368      numberNonZero++;
     369      int otherRow = parent_[iPivot];
     370      if (sign_[iPivot]<0) {
     371        array[iPivot] = - pivotValue;
     372        array[otherRow] += pivotValue;
     373      } else {
     374        array[otherRow] += pivotValue;
     375      }
     376    }
     377    iPivot = toRoot_[iPivot];
     378  }
     379  array[numberRows_]=0.0;
     380  return numberNonZero;
    66381}
    67382/* Updates one column (FTRAN) to/from array
    68    number returned is negative if no room
    69    ** For large problems you should ALWAYS know where the nonzeros
     383    ** For large problems you should ALWAYS know where the nonzeros
    70384   are, so please try and migrate to previous method after you
    71385   have got code working using this simple method - thank you!
    72    (the only exception is if you know input is dense e.g. rhs)
    73    region starts as zero and is zero at end */
     386   (the only exception is if you know input is dense e.g. rhs) */
    74387int
    75 ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse,
    76                         double array[] ) const
    77 {
    78   abort();
    79   return 1;
     388ClpNetworkBasis::updateColumn ( double array[] ) const
     389{
     390  int iPivot=leaf_;
     391  int numberNonZero=0;
     392  while (iPivot!=numberRows_) {
     393    double pivotValue = array[iPivot];
     394    if (pivotValue) {
     395      numberNonZero++;
     396      int otherRow = parent_[iPivot];
     397      if (sign_[iPivot]<0) {
     398        array[iPivot] = - pivotValue;
     399        array[otherRow] += pivotValue;
     400      } else {
     401        array[otherRow] += pivotValue;
     402      }
     403    }
     404    iPivot = toRoot_[iPivot];
     405  }
     406  array[numberRows_]=0.0;
     407  return numberNonZero;
    80408}
    81409/* Updates one column transpose (BTRAN)
     
    86414   returns number of nonzeros */
    87415int
    88 ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse,
    89                                           double array[] ) const
     416ClpNetworkBasis::updateColumnTranspose ( double array[] ) const
    90417{
    91418  abort();
    92419  return 1;
    93420}
    94 /* Updates one column (BTRAN) from region2
    95    region1 starts as zero and is zero at end */
     421/* Updates one column (BTRAN) from region2 */
    96422int
    97 ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse,
    98                           CoinIndexedVector * regionSparse2) const
     423ClpNetworkBasis::updateColumnTranspose (
     424                                        CoinIndexedVector * regionSparse2) const
    99425{
    100426  abort();
  • branches/unlabeled-1.1.2/ClpNetworkMatrix.cpp

    r118 r119  
    544544      be modified if doing column generation */
    545545void
    546 ClpNetworkMatrix::unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
     546ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
    547547                   int iColumn) const
    548548{
  • branches/unlabeled-1.1.2/ClpPlusMinusOneMatrix.cpp

    r118 r119  
    430430    }
    431431  } else if (numberInRowArray==2) {
    432     // do by rows when two rows
    433     int iRow;
     432    // do by rows when two rows (do longer first)
     433    int iRow0 = whichRow[0];
     434    int iRow1 = whichRow[1];
     435    if (startPositive[iRow0+1]-startPositive[iRow0]<
     436        startPositive[iRow1+1]-startPositive[iRow1]) {
     437      int temp = iRow0;
     438      iRow0 = iRow1;
     439      iRow1 = temp;
     440    }
    434441    int numberOriginal;
    435442    int i;
    436443    numberNonZero=0;
    437444    double value;
    438     iRow = whichRow[0];
    439     value = pi[iRow]*scalar;
     445    value = pi[iRow0]*scalar;
    440446    CoinBigIndex j;
    441     for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     447    for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
    442448      int iColumn = column[j];
    443449      index[numberNonZero++]=iColumn;
    444450      array[iColumn] = value;
    445451    }
    446     for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     452    for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
    447453      int iColumn = column[j];
    448454      index[numberNonZero++]=iColumn;
    449455      array[iColumn] = -value;
    450456    }
    451     iRow = whichRow[1];
    452     value = pi[iRow]*scalar;
    453     for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     457    value = pi[iRow1]*scalar;
     458    for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
    454459      int iColumn = column[j];
    455460      double value2= array[iColumn];
     
    462467      array[iColumn] = value2;
    463468    }
    464     for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     469    for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
    465470      int iColumn = column[j];
    466471      double value2= array[iColumn];
     
    595600      be modified if doing column generation */
    596601void
    597 ClpPlusMinusOneMatrix::unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
    598                    int iColumn) const
     602ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
     603                              CoinIndexedVector * rowArray,
     604                              int iColumn) const
    599605{
    600606  CoinBigIndex j=startPositive_[iColumn];
  • branches/unlabeled-1.1.2/include/ClpNetworkBasis.hpp

    r118 r119  
    1313class ClpMatrixBase;
    1414class CoinIndexedVector;
     15class ClpSimplex;
     16
    1517/** This deals with Factorization and Updates for network structures
    1618 */
     
    2527  /// Default constructor
    2628    ClpNetworkBasis (  );
     29  /// Constructor from CoinFactorization
     30  ClpNetworkBasis(const ClpSimplex * model,
     31                  int numberRows, const double * pivotRegion,
     32                  const int * permuteBack,const int * startColumn,
     33                  const int * numberInColumn,
     34                  const int * indexRow, const double * element);
    2735  /// Copy constructor
    2836  ClpNetworkBasis ( const ClpNetworkBasis &other);
     
    6169   which user may want to know about */
    6270  //@{
    63   /** Updates one column (FTRAN) from region2
    64       number returned is negative if no room
    65       region1 starts as zero and is zero at end */
    66   int updateColumn ( CoinIndexedVector * regionSparse,
    67                      CoinIndexedVector * regionSparse2);
     71  /** Updates one column (FTRAN) from region */
     72  int updateColumn ( CoinIndexedVector * regionSparse2);
    6873  /** Updates one column (FTRAN) to/from array
    69       number returned is negative if no room
    7074      ** For large problems you should ALWAYS know where the nonzeros
    7175      are, so please try and migrate to previous method after you
    7276      have got code working using this simple method - thank you!
    73       (the only exception is if you know input is dense e.g. rhs)
    74       region starts as zero and is zero at end */
    75   int updateColumn ( CoinIndexedVector * regionSparse,
    76                         double array[] ) const;
     77      (the only exception is if you know input is dense e.g. rhs) */
     78  int updateColumn ( double array[] ) const;
    7779  /** Updates one column transpose (BTRAN)
    7880      ** For large problems you should ALWAYS know where the nonzeros
     
    8183      (the only exception is if you know input is dense e.g. dense objective)
    8284      returns number of nonzeros */
    83   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    84                                  double array[] ) const;
    85   /** Updates one column (BTRAN) from region2
    86       region1 starts as zero and is zero at end */
    87   int updateColumnTranspose ( CoinIndexedVector * regionSparse,
    88                               CoinIndexedVector * regionSparse2) const;
     85  int updateColumnTranspose ( double array[] ) const;
     86  /** Updates one column (BTRAN) from region2 */
     87  int updateColumnTranspose ( CoinIndexedVector * regionSparse2) const;
    8988  //@}
    9089////////////////// data //////////////////
     
    9998  /// Number of Columns in factorization
    10099  int numberColumns_;
    101   /// Maximum number of pivots before factorization
    102   int maximumPivots_;
    103   /// Number pivots since last factorization
    104   int numberPivots_;
    105   /// Pivot order for each Column
    106   int *pivotColumn_;
    107   /// Permutation vector for pivot row order
    108   int *permute_;
    109   /// DePermutation vector for pivot row order
    110   int *permuteBack_;
     100  /// Index of root
     101  int root_;
     102  /// Index of extreme leaf
     103  int leaf_;
     104  /// model
     105  const ClpSimplex * model_;
     106  /// Parent for each column
     107  int * parent_;
     108  /// Descendant
     109  int * descendant_;
     110  /// Pivot row
     111  int * pivot_;
     112  /// Right sibling
     113  int * rightSibling_;
     114  /// Left sibling
     115  int * leftSibling_;
     116  /// Sign of pivot
     117  double * sign_;
     118  /// Stack
     119  int * stack_;
     120  /// Next one towards leaf
     121  int * toLeaf_;
     122  /// Next one towards root
     123  int * toRoot_;
     124  /// To mark rows
     125  char * mark_;
    111126  //@}
    112127};
  • branches/unlabeled-1.1.2/include/ClpNetworkMatrix.hpp

    r118 r119  
    6565      Note that model is NOT const.  Bounds and objective could
    6666      be modified if doing column generation */
    67   virtual void unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
     67  virtual void unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
    6868                   int column) const ;
    6969  /** Adds multiple of a column into an CoinIndexedvector
  • branches/unlabeled-1.1.2/include/ClpPlusMinusOneMatrix.hpp

    r118 r119  
    6363      Note that model is NOT const.  Bounds and objective could
    6464      be modified if doing column generation */
    65   virtual void unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
     65  virtual void unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
    6666                   int column) const ;
    6767  /** Adds multiple of a column into an CoinIndexedvector
  • trunk/ClpFactorization.cpp

    r118 r119  
    55#include "ClpFactorization.hpp"
    66#include "CoinHelperFunctions.hpp"
     7#include "CoinIndexedVector.hpp"
    78#include "ClpSimplex.hpp"
    89#include "ClpMatrixBase.hpp"
    910#include "ClpNetworkBasis.hpp"
     11#include "ClpNetworkMatrix.hpp"
    1012
    1113
     
    7274                              double areaFactor )
    7375{
    74   // maybe for speed will be better to leave as many regions as possible
    75   gutsOfDestructor();
    76   gutsOfInitialize(2);
    77   if (areaFactor)
    78     areaFactor_ = areaFactor;
    79   int numberBasic = 0;
    80   CoinBigIndex numberElements=0;
    81   int numberRowBasic=0;
    82 
    83   // compute how much in basis
    84 
    85   int i;
    86 
    87   for (i=0;i<numberRows;i++) {
    88     if (rowIsBasic[i]>=0)
    89       numberRowBasic++;
    90   }
    91 
    92   numberBasic = numberRowBasic;
    93   for (i=0;i<numberColumns;i++) {
    94     if (columnIsBasic[i]>=0) {
    95       numberBasic++;
    96     }
    97   }
    98   numberElements += matrix->numberInBasis(columnIsBasic);
    99   if ( numberBasic > numberRows ) {
     76  if (!networkBasis_||1) {
     77    // see if matrix a network
     78    ClpNetworkMatrix* networkMatrix =
     79      dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
     80    // If network - still allow ordinary factorization first time for laziness
     81
     82    int saveMaximumPivots = maximumPivots();
     83    if (networkMatrix)
     84      maximumPivots(1);
     85    // maybe for speed will be better to leave as many regions as possible
     86    gutsOfDestructor();
     87    gutsOfInitialize(2);
     88    if (areaFactor)
     89      areaFactor_ = areaFactor;
     90    int numberBasic = 0;
     91    CoinBigIndex numberElements=0;
     92    int numberRowBasic=0;
     93   
     94    // compute how much in basis
     95   
     96    int i;
     97   
     98    for (i=0;i<numberRows;i++) {
     99      if (rowIsBasic[i]>=0)
     100        numberRowBasic++;
     101    }
     102   
     103    numberBasic = numberRowBasic;
     104    for (i=0;i<numberColumns;i++) {
     105      if (columnIsBasic[i]>=0) {
     106        numberBasic++;
     107      }
     108    }
     109    numberElements += matrix->numberInBasis(columnIsBasic);
     110    if ( numberBasic > numberRows ) {
    100111    return -2; // say too many in basis
    101   }
    102   numberElements = 3 * numberBasic + 3 * numberElements + 10000;
    103   getAreas ( numberRows, numberBasic, numberElements,
    104              2 * numberElements );
    105   //fill
    106   //copy
    107   numberBasic=0;
    108   numberElements=0;
    109   for (i=0;i<numberRows;i++) {
    110     if (rowIsBasic[i]>=0) {
    111       indexRowU_[numberElements]=i;
    112       indexColumnU_[numberElements]=numberBasic;
    113       elementU_[numberElements++]=slackValue_;
    114       numberBasic++;
    115     }
    116   }
    117   numberElements +=matrix->fillBasis(model, columnIsBasic, numberBasic,
    118                                      indexRowU_+numberElements,
    119                                      indexColumnU_+numberElements,
    120                                      elementU_+numberElements);
    121   lengthU_ = numberElements;
    122 
    123   preProcess ( 0 );
    124   factor (  );
    125   numberBasic=0;
    126   if (status_ == 0) {
    127     int * permuteBack = permuteBack_;
    128     int * back = pivotColumnBack_;
     112    }
     113    numberElements = 3 * numberBasic + 3 * numberElements + 10000;
     114    getAreas ( numberRows, numberBasic, numberElements,
     115               2 * numberElements );
     116    //fill
     117    //copy
     118    numberBasic=0;
     119    numberElements=0;
    129120    for (i=0;i<numberRows;i++) {
    130121      if (rowIsBasic[i]>=0) {
    131         rowIsBasic[i]=permuteBack[back[numberBasic++]];
    132       }
    133     }
    134     for (i=0;i<numberColumns;i++) {
    135       if (columnIsBasic[i]>=0) {
    136         columnIsBasic[i]=permuteBack[back[numberBasic++]];
    137       }
    138     }
    139     if (increasingRows_>1) {
    140       // Set up permutation vector
    141       if (increasingRows_<3) {
    142         // these arrays start off as copies of permute
    143         // (and we could use permute_ instead of pivotColumn (not back though))
    144         ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
    145         ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
    146       }
    147     } else {
    148       // Set up permutation vector
    149       // (we could use permute_ instead of pivotColumn (not back though))
     122        indexRowU_[numberElements]=i;
     123        indexColumnU_[numberElements]=numberBasic;
     124        elementU_[numberElements++]=slackValue_;
     125        numberBasic++;
     126      }
     127    }
     128    numberElements +=matrix->fillBasis(model, columnIsBasic, numberBasic,
     129                                       indexRowU_+numberElements,
     130                                       indexColumnU_+numberElements,
     131                                       elementU_+numberElements);
     132    lengthU_ = numberElements;
     133   
     134    preProcess ( 0 );
     135    factor (  );
     136    numberBasic=0;
     137    if (status_ == 0) {
     138      int * permuteBack = permuteBack_;
     139      int * back = pivotColumnBack_;
     140      for (i=0;i<numberRows;i++) {
     141        if (rowIsBasic[i]>=0) {
     142          rowIsBasic[i]=permuteBack[back[numberBasic++]];
     143        }
     144      }
     145      for (i=0;i<numberColumns;i++) {
     146        if (columnIsBasic[i]>=0) {
     147          columnIsBasic[i]=permuteBack[back[numberBasic++]];
     148        }
     149      }
     150      if (increasingRows_>1) {
     151        // Set up permutation vector
     152        if (increasingRows_<3) {
     153          // these arrays start off as copies of permute
     154          // (and we could use permute_ instead of pivotColumn (not back though))
     155          ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
     156          ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
     157        }
     158      } else {
     159        // Set up permutation vector
     160        // (we could use permute_ instead of pivotColumn (not back though))
     161        for (i=0;i<numberRows_;i++) {
     162          int k=pivotColumn_[i];
     163          pivotColumn_[i]=pivotColumnBack_[i];
     164          pivotColumnBack_[i]=k;
     165        }
     166      }
     167    } else if (status_ == -1) {
     168      // mark as basic or non basic
    150169      for (i=0;i<numberRows_;i++) {
    151         int k=pivotColumn_[i];
    152         pivotColumn_[i]=pivotColumnBack_[i];
    153         pivotColumnBack_[i]=k;
    154       }
    155     }
    156   } else if (status_ == -1) {
    157     // mark as basic or non basic
    158     for (i=0;i<numberRows_;i++) {
    159       if (rowIsBasic[i]>=0) {
    160         if (pivotColumn_[numberBasic]>=0)
    161           rowIsBasic[i]=pivotColumn_[numberBasic];
    162         else
    163           rowIsBasic[i]=-1;
    164         numberBasic++;
    165       }
    166     }
    167     for (i=0;i<numberColumns;i++) {
    168       if (columnIsBasic[i]>=0) {
    169         if (pivotColumn_[numberBasic]>=0)
    170           columnIsBasic[i]=pivotColumn_[numberBasic];
    171         else
    172           columnIsBasic[i]=-1;
    173         numberBasic++;
    174       }
    175     }
     170        if (rowIsBasic[i]>=0) {
     171          if (pivotColumn_[numberBasic]>=0)
     172            rowIsBasic[i]=pivotColumn_[numberBasic];
     173          else
     174            rowIsBasic[i]=-1;
     175          numberBasic++;
     176        }
     177      }
     178      for (i=0;i<numberColumns;i++) {
     179        if (columnIsBasic[i]>=0) {
     180          if (pivotColumn_[numberBasic]>=0)
     181            columnIsBasic[i]=pivotColumn_[numberBasic];
     182          else
     183            columnIsBasic[i]=-1;
     184          numberBasic++;
     185        }
     186      }
     187    }
     188    if (networkMatrix) {
     189      maximumPivots(saveMaximumPivots);
     190      if (!status_) {
     191        // create network factorization
     192        delete networkBasis_; // temp
     193        networkBasis_ = new ClpNetworkBasis(model,numberRows_,
     194                                            pivotRegion_,
     195                                            permuteBack_,
     196                                            startColumnU_,
     197                                            numberInColumn_,
     198                                            indexRowU_,
     199                                            elementU_);
     200        // kill off arrays in ordinary factorization
     201        //gutsOfDestructor(); temp
     202      }
     203    }
     204  } else {
     205    // network - fake factorization - do nothing
     206    status_=0;
    176207  }
    177208
     
    197228                                            checkBeforeModifying);
    198229  } else {
    199     return networkBasis_->replaceColumn(regionSparse,
     230    int returnCode = CoinFactorization::replaceColumn(regionSparse,
     231                                            pivotRow,
     232                                            pivotCheck,
     233                                            checkBeforeModifying);
     234    networkBasis_->replaceColumn(regionSparse,
    200235                                        pivotRow);
     236    return returnCode;
    201237  }
    202238}
     
    215251                                           FTUpdate);
    216252  } else {
    217     return networkBasis_->updateColumn(regionSparse,
    218                                        regionSparse2);
     253    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
     254    int returnCode = CoinFactorization::updateColumn(regionSparse,
     255                                           regionSparse2, FTUpdate);
     256    networkBasis_->updateColumn(save);
     257    int i;
     258    double * array = regionSparse2->denseVector();
     259    double * array2 = save->denseVector();
     260    for (i=0;i<numberRows_;i++)
     261      assert (array2[i]==array[i]);
     262    delete save;
     263    return returnCode;
    219264  }
    220265}
     
    234279                                           array);
    235280  } else {
    236     return networkBasis_->updateColumn(regionSparse,
    237                                                     array);
     281    double * save = new double [numberRows_+1];
     282    memcpy(save,array,(numberRows_+1)*sizeof(double));
     283    int returnCode = CoinFactorization::updateColumn(regionSparse,
     284                                           array);
     285    networkBasis_->updateColumn(save);
     286    int i;
     287    for (i=0;i<numberRows_;i++)
     288      assert (save[i]==array[i]);
     289    delete [] save;
     290    return returnCode;
    238291  }
    239292}
     
    248301                                          double array[] ) const
    249302{
    250   if (!networkBasis_) {
     303  if (!networkBasis_||1) {
    251304    return CoinFactorization::updateColumnTranspose(regionSparse,
    252305                                                    array);
    253306  } else {
    254     return networkBasis_->updateColumnTranspose(regionSparse,
    255                                                     array);
     307    return networkBasis_->updateColumnTranspose(array);
    256308  }
    257309}
     
    262314                          CoinIndexedVector * regionSparse2) const
    263315{
    264   if (!networkBasis_) {
     316  if (!networkBasis_||1) {
    265317    return CoinFactorization::updateColumnTranspose(regionSparse,
    266318                                                    regionSparse2);
    267319  } else {
    268     return networkBasis_->updateColumnTranspose(regionSparse,
    269                                                     regionSparse2);
    270   }
    271 }
     320    return networkBasis_->updateColumnTranspose( regionSparse2);
     321  }
     322}
     323/* makes a row copy of L for speed and to allow very sparse problems */
     324void
     325ClpFactorization::goSparse()
     326{
     327  if (!networkBasis_)
     328    CoinFactorization::goSparse();
     329}
  • trunk/ClpModel.cpp

    r114 r119  
    2424#include "CoinMpsIO.hpp"
    2525#include "ClpMessage.hpp"
     26#include "ClpLinearObjective.hpp"
    2627
    2728static double cpuTime()
     
    117118  delete [] columnLower_;
    118119  delete [] columnUpper_;
    119   delete [] objective_;
     120  delete objective_;
    120121  columnLower_=NULL;
    121122  columnUpper_=NULL;
     
    169170  rowLower_=ClpCopyOfArray(rowlb,numberRows_,-DBL_MAX);
    170171  rowUpper_=ClpCopyOfArray(rowub,numberRows_,DBL_MAX);
    171   objective_=ClpCopyOfArray(obj,numberColumns_,0.0);
     172  double * objective=ClpCopyOfArray(obj,numberColumns_,0.0);
     173  objective_ = new ClpLinearObjective(objective,numberColumns_);
     174  delete [] objective;
    172175  rowObjective_=ClpCopyOfArray(rowObjective,numberRows_);
    173176  columnLower_=ClpCopyOfArray(collb,numberColumns_,0.0);
     
    362365    columnLower_ = ClpCopyOfArray ( rhs.columnLower_, numberColumns_ );
    363366    columnUpper_ = ClpCopyOfArray ( rhs.columnUpper_, numberColumns_ );
    364     objective_ = ClpCopyOfArray ( rhs.objective_, numberColumns_ );
     367    if (rhs.objective_)
     368      objective_  = rhs.objective_->clone();
     369    else
     370      objective_ = NULL;
    365371    rowObjective_ = ClpCopyOfArray ( rhs.rowObjective_, numberRows_ );
    366372    status_ = ClpCopyOfArray( rhs.status_,numberColumns_+numberRows_);
     
    627633  reducedCost_ = resizeDouble(reducedCost_,numberColumns_,
    628634                              newNumberColumns,0.0,true);
    629   objective_ = resizeDouble(objective_,numberColumns_,
    630                             newNumberColumns,0.0,true);
     635  objective_->resize(newNumberColumns);
    631636  columnLower_ = resizeDouble(columnLower_,numberColumns_,
    632637                              newNumberColumns,0.0,true);
     
    726731  reducedCost_ = deleteDouble(reducedCost_,numberColumns_,
    727732                              number, which, newSize);
    728   objective_ = deleteDouble(objective_,numberColumns_,
    729                               number, which, newSize);
     733  objective_->deleteSome(number, which);
    730734  columnLower_ = deleteDouble(columnLower_,numberColumns_,
    731735                              number, which, newSize);
  • trunk/ClpPackedMatrix.cpp

    r118 r119  
    730730      be modified if doing column generation */
    731731void
    732 ClpPackedMatrix::unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
     732ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
    733733                   int iColumn) const
    734734{
  • trunk/ClpSimplex.cpp

    r115 r119  
    2020#include "ClpNonLinearCost.hpp"
    2121#include "ClpMessage.hpp"
     22#include "ClpLinearObjective.hpp"
    2223#include <cfloat>
    2324
     
    243244  CoinIndexedVector  * workSpace = rowArray_[0];
    244245
    245   double * array = new double [numberRows_];
     246  double * array = new double [numberRows_+1]; // +1 for network
    246247  double * save = new double [numberRows_];
    247248  double * previous = new double [numberRows_];
     
    343344  CoinIndexedVector  * workSpace = rowArray_[0];
    344345
    345   double * array = new double [numberRows_];
     346  double * array = new double [numberRows_+1]; // +1 for network
    346347  double * save = new double [numberRows_];
    347348  double * previous = new double [numberRows_];
     
    16411642*/
    16421643void
    1643 ClpSimplex::unpack(CoinIndexedVector * rowArray)
     1644ClpSimplex::unpack(CoinIndexedVector * rowArray) const
    16441645{
    16451646  rowArray->clear();
     
    16531654}
    16541655void
    1655 ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence)
     1656ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence) const
    16561657{
    16571658  rowArray->clear();
     
    17261727    objectiveWork_ = cost_;
    17271728    rowObjectiveWork_ = cost_+numberColumns_;
    1728     memcpy(objectiveWork_,objective_,numberColumns_*sizeof(double));
     1729    memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double));
    17291730    if (rowObjective_)
    17301731      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
     
    24282429    if (outDoubleArray(rowUpper_,numberRows_,fp))
    24292430        return 1;
    2430     if (outDoubleArray(objective_,numberColumns_,fp))
     2431    if (outDoubleArray(objective(),numberColumns_,fp))
    24312432        return 1;
    24322433    if (outDoubleArray(rowObjective_,numberRows_,fp))
     
    26162617    if (inDoubleArray(rowUpper_,numberRows_,fp))
    26172618        return 1;
    2618     if (inDoubleArray(objective_,numberColumns_,fp))
     2619    double * objective;
     2620    if (inDoubleArray(objective,numberColumns_,fp))
    26192621        return 1;
     2622    delete objective_;
     2623    objective_ = new ClpLinearObjective(objective,numberColumns_);
     2624    delete [] objective;
    26202625    if (inDoubleArray(rowObjective_,numberRows_,fp))
    26212626        return 1;
     
    32303235    double * dj = new double [numberColumns_];
    32313236    double * solution = columnActivity_;
     3237    const double * linearObjective = objective();
    32323238    //double objectiveValue=0.0;
    32333239    int iColumn;
    32343240    for (iColumn=0;iColumn<numberColumns_;iColumn++)
    3235       dj[iColumn] = optimizationDirection_*objective_[iColumn];
     3241      dj[iColumn] = optimizationDirection_*linearObjective[iColumn];
    32363242    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    32373243      // assume natural place is closest to zero
  • trunk/ClpSimplexDual.cpp

    r111 r119  
    1 
     1 
    22// Copyright (C) 2002, International Business Machines
    33// Corporation and others.  All Rights Reserved.
     
    100100#include "CoinWarmStartBasis.hpp"
    101101#include "ClpDualRowDantzig.hpp"
     102#include "ClpPlusMinusOneMatrix.hpp"
    102103#include "ClpMessage.hpp"
    103104#include <cfloat>
     
    196197
    197198  */
    198 
    199199  algorithm_ = -1;
    200200
     
    18621862  if (!rowScale_&&(handler_->logLevel()&32)) {
    18631863    double * objectiveSimplex
    1864       = ClpCopyOfArray(objective_,numberColumns_,0.0);
     1864      = ClpCopyOfArray(objective(),numberColumns_,0.0);
    18651865    double * rowObjectiveSimplex
    18661866      = ClpCopyOfArray(rowObjective_,numberRows_,0.0);
  • trunk/Makefile.Clp

    r118 r119  
    2121LIBSRC += ClpNetworkMatrix.cpp
    2222LIBSRC += ClpNonLinearCost.cpp
     23LIBSRC += ClpObjective.cpp
     24LIBSRC += ClpLinearObjective.cpp
    2325LIBSRC += ClpPackedMatrix.cpp
    2426LIBSRC += ClpPlusMinusOneMatrix.cpp
  • trunk/Samples/Makefile.network

    r68 r119  
    5959        LDFLAGS += -lhistory -lreadline -ltermcap
    6060endif
     61LDFLAGS += -lefence
    6162###############################################################################
    6263
  • trunk/Samples/network.cpp

    r68 r119  
    1212#include "ClpSimplex.hpp"
    1313#include "ClpFactorization.hpp"
     14#include "ClpNetworkMatrix.hpp"
    1415#include <stdio.h>
    1516#include <assert.h>
     
    139140    objective[i]=cost[i];
    140141  }
     142    // Create Packed Matrix
     143  CoinPackedMatrix matrix;
     144  int * lengths = NULL;
     145  matrix.assignMatrix(true,numberRows,numberColumns,
     146                      2*numberColumns,element,row,start,lengths);
     147  ClpNetworkMatrix network(matrix);
    141148  // load model
    142149
    143   model.loadProblem(numberColumns,numberRows,start,row,element,
     150  model.loadProblem(network,
    144151                    lowerColumn,upperColumn,objective,
    145152                    lower,upper);
     
    158165  delete [] row;
    159166  model.factorization()->maximumPivots(200+model.numberRows()/100);
     167  model.factorization()->maximumPivots(2);
     168  model.messageHandler()->setLogLevel(63);
    160169  model.dual();
    161170  return 0;
  • trunk/include/ClpFactorization.hpp

    r118 r119  
    106106  //@}
    107107   
     108  /**@name other stuff */
     109  //@{
     110  /** makes a row copy of L for speed and to allow very sparse problems */
     111  void goSparse();
     112  //@}
     113
    108114////////////////// data //////////////////
    109115private:
  • trunk/include/ClpMatrixBase.hpp

    r118 r119  
    8989      Note that model is NOT const.  Bounds and objective could
    9090      be modified if doing column generation (just for this variable) */
    91   virtual void unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
     91  virtual void unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
    9292                   int column) const =0;
    9393  /** Purely for column generation and similar ideas.  Allows
  • trunk/include/ClpModel.hpp

    r114 r119  
    1313#include "CoinMessageHandler.hpp"
    1414#include "ClpParameters.hpp"
     15#include "ClpObjective.hpp"
    1516
    1617#define CLP_INFINITY 1e30
     
    205206   inline const double* getRowUpper() const     { return rowUpper_; }
    206207   /// Objective
    207    inline double * objective() const            { return objective_; }
    208    inline const double * getObjCoefficients() const { return objective_; }
     208   inline double * objective() const           
     209  { double offset; return objective_->gradient(NULL,offset); }
     210   inline const double * getObjCoefficients() const
     211  { double offset; return objective_->gradient(NULL,offset); }
    209212   /// Row Objective
    210213   inline double * rowObjective() const         { return rowObjective_; }
     
    399402  double* rowUpper_;
    400403  /// Objective
    401   double * objective_;
     404  ClpObjective * objective_;
    402405  /// Row Objective (? sign)  - may be NULL
    403406  double * rowObjective_;
  • trunk/include/ClpPackedMatrix.hpp

    r118 r119  
    8989      Note that model is NOT const.  Bounds and objective could
    9090      be modified if doing column generation */
    91   virtual void unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
     91  virtual void unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
    9292                   int column) const ;
    9393  /** Adds multiple of a column into an CoinIndexedvector
  • trunk/include/ClpSimplex.hpp

    r115 r119  
    330330     Also applies scaling if needed
    331331  */
    332   void unpack(CoinIndexedVector * rowArray);
     332  void unpack(CoinIndexedVector * rowArray) const ;
    333333  /**
    334334     Unpacks one column of the matrix into indexed array
     
    336336     Also applies scaling if needed
    337337  */
    338   void unpack(CoinIndexedVector * rowArray,int sequence);
     338  void unpack(CoinIndexedVector * rowArray,int sequence) const;
    339339 
    340340  /**
Note: See TracChangeset for help on using the changeset viewer.