source: trunk/ClpFactorization.cpp @ 118

Last change on this file since 118 was 118, checked in by forrest, 17 years ago

Adding Network matrix and PlusMinusOne?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include "ClpFactorization.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "ClpSimplex.hpp"
8#include "ClpMatrixBase.hpp"
9#include "ClpNetworkBasis.hpp"
10
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpFactorization::ClpFactorization () :
20   CoinFactorization() 
21{
22  networkBasis_ = NULL;
23}
24
25//-------------------------------------------------------------------
26// Copy constructor
27//-------------------------------------------------------------------
28ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
29   CoinFactorization(rhs) 
30{
31  if (rhs.networkBasis_)
32    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
33  else
34    networkBasis_=NULL;
35}
36
37ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
38   CoinFactorization(rhs) 
39{
40  networkBasis_=NULL;
41}
42
43//-------------------------------------------------------------------
44// Destructor
45//-------------------------------------------------------------------
46ClpFactorization::~ClpFactorization () 
47{
48  delete networkBasis_;
49}
50
51//----------------------------------------------------------------
52// Assignment operator
53//-------------------------------------------------------------------
54ClpFactorization &
55ClpFactorization::operator=(const ClpFactorization& rhs)
56{
57  if (this != &rhs) {
58    CoinFactorization::operator=(rhs);
59    delete networkBasis_;
60    if (rhs.networkBasis_)
61      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
62    else
63      networkBasis_=NULL;
64  }
65  return *this;
66}
67int 
68ClpFactorization::factorize ( const ClpSimplex * model,
69                              const ClpMatrixBase * matrix, 
70                              int numberRows, int numberColumns,
71                              int rowIsBasic[], int columnIsBasic[] , 
72                              double areaFactor )
73{
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 ) {
100    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_;
129    for (i=0;i<numberRows;i++) {
130      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))
150      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    }
176  }
177
178  return status_;
179}
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}
Note: See TracBrowser for help on using the repository browser.