source: trunk/ClpFactorization.cpp @ 138

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

Fix for odd no btran case

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.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 "CoinIndexedVector.hpp"
8#include "ClpSimplex.hpp"
9#include "ClpMatrixBase.hpp"
10#include "ClpNetworkBasis.hpp"
11#include "ClpNetworkMatrix.hpp"
12//#define CHECK_NETWORK
13#ifdef CHECK_NETWORK
14const static bool doCheck=true;
15#else
16const static bool doCheck=false;
17#endif
18
19//#############################################################################
20// Constructors / Destructor / Assignment
21//#############################################################################
22
23//-------------------------------------------------------------------
24// Default Constructor
25//-------------------------------------------------------------------
26ClpFactorization::ClpFactorization () :
27   CoinFactorization() 
28{
29  networkBasis_ = NULL;
30}
31
32//-------------------------------------------------------------------
33// Copy constructor
34//-------------------------------------------------------------------
35ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
36   CoinFactorization(rhs) 
37{
38  if (rhs.networkBasis_)
39    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
40  else
41    networkBasis_=NULL;
42}
43
44ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
45   CoinFactorization(rhs) 
46{
47  networkBasis_=NULL;
48}
49
50//-------------------------------------------------------------------
51// Destructor
52//-------------------------------------------------------------------
53ClpFactorization::~ClpFactorization () 
54{
55  delete networkBasis_;
56}
57
58//----------------------------------------------------------------
59// Assignment operator
60//-------------------------------------------------------------------
61ClpFactorization &
62ClpFactorization::operator=(const ClpFactorization& rhs)
63{
64  if (this != &rhs) {
65    CoinFactorization::operator=(rhs);
66    delete networkBasis_;
67    if (rhs.networkBasis_)
68      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
69    else
70      networkBasis_=NULL;
71  }
72  return *this;
73}
74int 
75ClpFactorization::factorize ( const ClpSimplex * model,
76                              const ClpMatrixBase * matrix, 
77                              int numberRows, int numberColumns,
78                              int rowIsBasic[], int columnIsBasic[] , 
79                              double areaFactor )
80{
81  if (!networkBasis_||doCheck) {
82    // see if matrix a network
83    ClpNetworkMatrix* networkMatrix =
84      dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
85    // If network - still allow ordinary factorization first time for laziness
86    int saveMaximumPivots = maximumPivots();
87    delete networkBasis_;
88    networkBasis_ = NULL;
89    if (networkMatrix&&!doCheck)
90      maximumPivots(1);
91    // maybe for speed will be better to leave as many regions as possible
92    gutsOfDestructor();
93    gutsOfInitialize(2);
94    if (areaFactor)
95      areaFactor_ = areaFactor;
96    int numberBasic = 0;
97    CoinBigIndex numberElements=0;
98    int numberRowBasic=0;
99   
100    // compute how much in basis
101   
102    int i;
103   
104    for (i=0;i<numberRows;i++) {
105      if (rowIsBasic[i]>=0)
106        numberRowBasic++;
107    }
108   
109    numberBasic = numberRowBasic;
110    for (i=0;i<numberColumns;i++) {
111      if (columnIsBasic[i]>=0) {
112        numberBasic++;
113      }
114    }
115    numberElements += matrix->numberInBasis(columnIsBasic);
116    if ( numberBasic > numberRows ) {
117    return -2; // say too many in basis
118    }
119    numberElements = 3 * numberBasic + 3 * numberElements + 10000;
120    getAreas ( numberRows, numberBasic, numberElements,
121               2 * numberElements );
122    //fill
123    //copy
124    numberBasic=0;
125    numberElements=0;
126    for (i=0;i<numberRows;i++) {
127      if (rowIsBasic[i]>=0) {
128        indexRowU_[numberElements]=i;
129        indexColumnU_[numberElements]=numberBasic;
130        elementU_[numberElements++]=slackValue_;
131        numberBasic++;
132      }
133    }
134    numberElements +=matrix->fillBasis(model, columnIsBasic, numberBasic, 
135                                       indexRowU_+numberElements, 
136                                       indexColumnU_+numberElements,
137                                       elementU_+numberElements);
138    lengthU_ = numberElements;
139   
140    preProcess ( 0 );
141    factor (  );
142    numberBasic=0;
143    if (status_ == 0) {
144      int * permuteBack = permuteBack_;
145      int * back = pivotColumnBack_;
146      for (i=0;i<numberRows;i++) {
147        if (rowIsBasic[i]>=0) {
148          rowIsBasic[i]=permuteBack[back[numberBasic++]];
149        }
150      }
151      for (i=0;i<numberColumns;i++) {
152        if (columnIsBasic[i]>=0) {
153          columnIsBasic[i]=permuteBack[back[numberBasic++]];
154        }
155      }
156      if (increasingRows_>1) {
157        // Set up permutation vector
158        if (increasingRows_<3) {
159          // these arrays start off as copies of permute
160          // (and we could use permute_ instead of pivotColumn (not back though))
161          ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
162          ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
163        }
164      } else {
165        // Set up permutation vector
166        // (we could use permute_ instead of pivotColumn (not back though))
167        for (i=0;i<numberRows_;i++) {
168          int k=pivotColumn_[i];
169          pivotColumn_[i]=pivotColumnBack_[i];
170          pivotColumnBack_[i]=k;
171        }
172      }
173    } else if (status_ == -1) {
174      // mark as basic or non basic
175      for (i=0;i<numberRows_;i++) {
176        if (rowIsBasic[i]>=0) {
177          if (pivotColumn_[numberBasic]>=0) 
178            rowIsBasic[i]=pivotColumn_[numberBasic];
179          else
180            rowIsBasic[i]=-1;
181          numberBasic++;
182        }
183      }
184      for (i=0;i<numberColumns;i++) {
185        if (columnIsBasic[i]>=0) {
186          if (pivotColumn_[numberBasic]>=0) 
187            columnIsBasic[i]=pivotColumn_[numberBasic];
188          else
189            columnIsBasic[i]=-1;
190          numberBasic++;
191        }
192      }
193    }
194    if (networkMatrix) {
195      maximumPivots(saveMaximumPivots);
196      if (!status_) {
197        // create network factorization
198        if (doCheck)
199          delete networkBasis_; // temp
200        networkBasis_ = new ClpNetworkBasis(model,numberRows_,
201                                            pivotRegion_,
202                                            permuteBack_,
203                                            startColumnU_,
204                                            numberInColumn_,
205                                            indexRowU_,
206                                            elementU_);
207        // kill off arrays in ordinary factorization
208        if (!doCheck)
209          gutsOfDestructor();
210      }
211    }
212    if (!status_&&!networkBasis_) {
213      // See if worth going sparse and when
214      if (numberFtranCounts_>100) {
215        ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
216        ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
217        ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
218        if (btranCountInput_) {
219          btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
220          btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
221          btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
222        } else {
223          // odd - we have not done any btrans
224          btranAverageAfterU_ = 1.0;
225          btranAverageAfterR_ = 1.0;
226          btranAverageAfterL_ = 1.0;
227        }
228      }
229      // scale back
230
231      ftranCountInput_ *= 0.8;
232      ftranCountAfterL_ *= 0.8;
233      ftranCountAfterR_ *= 0.8;
234      ftranCountAfterU_ *= 0.8;
235      btranCountInput_ *= 0.8;
236      btranCountAfterU_ *= 0.8;
237      btranCountAfterR_ *= 0.8;
238      btranCountAfterL_ *= 0.8;
239    }
240  } else {
241    // network - fake factorization - do nothing
242    status_=0;
243  }
244
245  return status_;
246}
247/* Replaces one Column to basis,
248   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
249   If checkBeforeModifying is true will do all accuracy checks
250   before modifying factorization.  Whether to set this depends on
251   speed considerations.  You could just do this on first iteration
252   after factorization and thereafter re-factorize
253   partial update already in U */
254int 
255ClpFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
256                      int pivotRow,
257                      double pivotCheck ,
258                      bool checkBeforeModifying)
259{
260  if (!networkBasis_) {
261    return CoinFactorization::replaceColumn(regionSparse,
262                                            pivotRow,
263                                            pivotCheck,
264                                            checkBeforeModifying);
265  } else {
266    if (doCheck) {
267      int returnCode = CoinFactorization::replaceColumn(regionSparse,
268                                                        pivotRow,
269                                                        pivotCheck,
270                                                        checkBeforeModifying);
271      networkBasis_->replaceColumn(regionSparse,
272                                   pivotRow);
273      return returnCode;
274    } else {
275      // increase number of pivots
276      numberPivots_++;
277      return networkBasis_->replaceColumn(regionSparse,
278                                   pivotRow);
279    }
280  }
281}
282
283/* Updates one column (FTRAN) from region2
284   number returned is negative if no room
285   region1 starts as zero and is zero at end */
286int 
287ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
288                                 CoinIndexedVector * regionSparse2,
289                                 bool FTUpdate) 
290{
291#ifdef CLP_DEBUG
292  regionSparse->checkClear();
293#endif
294  if (!networkBasis_) {
295    collectStatistics_ = true;
296    return CoinFactorization::updateColumn(regionSparse,
297                                           regionSparse2,
298                                           FTUpdate);
299    collectStatistics_ = false;
300  } else {
301#ifdef CHECK_NETWORK
302    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
303    int returnCode = CoinFactorization::updateColumn(regionSparse,
304                                                     regionSparse2, FTUpdate);
305    networkBasis_->updateColumn(regionSparse,save);
306    int i;
307    double * array = regionSparse2->denseVector();
308    double * array2 = save->denseVector();
309    for (i=0;i<numberRows_;i++) {
310      double value1 = array[i];
311      double value2 = array2[i];
312      assert (value1==value2);
313    }
314    delete save;
315    return returnCode;
316#else
317    return networkBasis_->updateColumn(regionSparse,regionSparse2);
318#endif
319  }
320}
321/* Updates one column (FTRAN) to/from array
322   number returned is negative if no room
323   ** For large problems you should ALWAYS know where the nonzeros
324   are, so please try and migrate to previous method after you
325   have got code working using this simple method - thank you!
326   (the only exception is if you know input is dense e.g. rhs)
327   region starts as zero and is zero at end */
328int 
329ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
330                        double array[] ) const
331{
332  if (!networkBasis_) {
333    return CoinFactorization::updateColumn(regionSparse,
334                                           array);
335  } else {
336#ifdef CHECK_NETWORK
337    double * save = new double [numberRows_+1];
338    memcpy(save,array,(numberRows_+1)*sizeof(double));
339    int returnCode = CoinFactorization::updateColumn(regionSparse,
340                                                     array);
341    networkBasis_->updateColumn(regionSparse, save);
342    int i;
343    for (i=0;i<numberRows_;i++)
344      assert (fabs(save[i]-array[i])<1.0e-8*(1.0+fabs(array[i])));
345    delete [] save;
346    return returnCode;
347#else
348    return networkBasis_->updateColumn(regionSparse, array);
349#endif
350  }
351}
352/* Updates one column transpose (BTRAN)
353   For large problems you should ALWAYS know where the nonzeros
354   are, so please try and migrate to previous method after you
355   have got code working using this simple method - thank you!
356   (the only exception is if you know input is dense e.g. dense objective)
357   returns number of nonzeros */
358int 
359ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
360                                          double array[] ) const
361{
362  if (!networkBasis_) {
363    return CoinFactorization::updateColumnTranspose(regionSparse,
364                                                    array);
365  } else {
366#ifdef CHECK_NETWORK
367    double * save = new double [numberRows_+1];
368    memcpy(save,array,(numberRows_+1)*sizeof(double));
369    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
370                                                              array);
371    networkBasis_->updateColumnTranspose(regionSparse, save);
372    int i;
373    for (i=0;i<numberRows_;i++)
374      assert (fabs(save[i]-array[i])<1.0e-8*(1.0+fabs(array[i])));
375    delete [] save;
376    return returnCode;
377#else
378    return networkBasis_->updateColumnTranspose(regionSparse, array);
379#endif
380  }
381}
382/* Updates one column (BTRAN) from region2
383   region1 starts as zero and is zero at end */
384int 
385ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
386                          CoinIndexedVector * regionSparse2) const
387{
388  if (!networkBasis_) {
389    collectStatistics_ = true;
390    return CoinFactorization::updateColumnTranspose(regionSparse,
391                                                    regionSparse2);
392    collectStatistics_ = false;
393  } else {
394#ifdef CHECK_NETWORK
395      CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
396      int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
397                                                                regionSparse2);
398      networkBasis_->updateColumnTranspose(regionSparse,save);
399      int i;
400      double * array = regionSparse2->denseVector();
401      double * array2 = save->denseVector();
402      for (i=0;i<numberRows_;i++) {
403        double value1 = array[i];
404        double value2 = array2[i];
405        assert (value1==value2);
406      }
407      delete save;
408      return returnCode;
409#else
410      return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
411#endif
412  }
413}
414/* makes a row copy of L for speed and to allow very sparse problems */
415void 
416ClpFactorization::goSparse()
417{
418  if (!networkBasis_) 
419    CoinFactorization::goSparse();
420}
421// Cleans up i.e. gets rid of network basis
422void 
423ClpFactorization::cleanUp()
424{
425  delete networkBasis_;
426  networkBasis_=NULL;
427  resetStatistics();
428}
429/// Says whether to redo pivot order
430bool 
431ClpFactorization::needToReorder() const
432{
433#ifdef CHECK_NETWORK
434  return true;
435#endif
436  if (!networkBasis_)
437    return true;
438  else
439    return false;
440}
Note: See TracBrowser for help on using the repository browser.