source: branches/pre/ClpFactorization.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago

new stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.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 ( ClpSimplex * model,
76                              int solveType, bool valuesPass)
77{
78  const ClpMatrixBase * matrix = model->clpMatrix(); 
79  int numberRows = model->numberRows();
80  int numberColumns = model->numberColumns();
81  // If too many compressions increase area
82  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
83    areaFactor_ *= 1.1;
84  }
85  if (!networkBasis_||doCheck) {
86    status_=-99;
87    int * pivotVariable = model->pivotVariable();
88    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
89    while (status_<-98) {
90     
91      int i;
92      int numberBasic=0;
93      int numberRowBasic;
94      for (i=0;i<numberRows;i++) {
95        if (model->getRowStatus(i) == ClpSimplex::basic) 
96          pivotVariable[numberBasic++]=i;
97      }
98      numberRowBasic=numberBasic;
99      for (i=0;i<numberColumns;i++) {
100        if (model->getColumnStatus(i) == ClpSimplex::basic) 
101          pivotVariable[numberBasic++]=i;
102      }
103      assert (numberBasic<=numberRows);
104      // see if matrix a network
105      ClpNetworkMatrix* networkMatrix =
106        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
107      // If network - still allow ordinary factorization first time for laziness
108      int saveMaximumPivots = maximumPivots();
109      delete networkBasis_;
110      networkBasis_ = NULL;
111      if (networkMatrix&&!doCheck)
112        maximumPivots(1);
113      while (status_==-99) {
114        // maybe for speed will be better to leave as many regions as possible
115        gutsOfDestructor();
116        gutsOfInitialize(2);
117        CoinBigIndex numberElements=numberRowBasic;
118
119        // compute how much in basis
120
121        int i;
122
123        numberElements +=matrix->fillBasis(model,
124                                           pivotVariable+numberRowBasic, 
125                                           numberRowBasic,
126                                           numberBasic-numberRowBasic,
127                                           NULL,NULL,NULL);
128
129        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
130        getAreas ( numberRows, numberBasic, numberElements,
131                   2 * numberElements );
132        //fill
133        //copy
134        numberElements=numberRowBasic;
135        for (i=0;i<numberRowBasic;i++) {
136          int iRow = pivotVariable[i];
137          indexRowU_[i]=iRow;
138          indexColumnU_[i]=i;
139          elementU_[i]=slackValue_;
140        }
141        numberElements +=matrix->fillBasis(model, 
142                                           pivotVariable+numberRowBasic, 
143                                           numberRowBasic, 
144                                           numberBasic-numberRowBasic,
145                                           indexRowU_+numberElements, 
146                                           indexColumnU_+numberElements,
147                                           elementU_+numberElements);
148        lengthU_ = numberElements;
149
150        preProcess ( 0 );
151        factor (  );
152        if (status_==-99) {
153          // get more memory
154          areaFactor(2.0*areaFactor());
155        }
156      }
157      // If we get here status is 0 or -1
158      if (status_ == 0) {
159        int * permuteBack = permuteBack_;
160        int * back = pivotColumnBack_;
161        int * pivotTemp = pivotColumn_;
162        ClpDisjointCopyN ( pivotVariable, numberRows_ , pivotTemp  );
163        // Redo pivot order
164        for (i=0;i<numberRowBasic;i++) {
165          int k = pivotTemp[i];
166          // so rowIsBasic[k] would be permuteBack[back[i]]
167          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
168        }
169        for (;i<numberRows;i++) {
170          int k = pivotTemp[i];
171          // so rowIsBasic[k] would be permuteBack[back[i]]
172          pivotVariable[permuteBack[back[i]]]=k;
173        }
174        // Set up permutation vector
175        // these arrays start off as copies of permute
176        // (and we could use permute_ instead of pivotColumn (not back though))
177        ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
178        ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
179        if (networkMatrix) {
180          maximumPivots(saveMaximumPivots);
181          // create network factorization
182          if (doCheck)
183            delete networkBasis_; // temp
184          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
185                                              pivotRegion_,
186                                              permuteBack_,
187                                              startColumnU_,
188                                              numberInColumn_,
189                                              indexRowU_,
190                                              elementU_);
191          // kill off arrays in ordinary factorization
192          if (!doCheck) {
193            gutsOfDestructor();
194#if 0
195            // but put back permute arrays so odd things will work
196            int numberRows = model->numberRows();
197            pivotColumnBack_ = new int [numberRows];
198            permute_ = new int [numberRows];
199            int i;
200            for (i=0;i<numberRows;i++) {
201              pivotColumnBack_[i]=i;
202              permute_[i]=i;
203            }
204#endif
205          }
206        } else {
207          // See if worth going sparse and when
208          if (numberFtranCounts_>100) {
209            ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
210            ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
211            ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
212            if (btranCountInput_) {
213              btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
214              btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
215              btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
216            } else {
217              // odd - we have not done any btrans
218              btranAverageAfterU_ = 1.0;
219              btranAverageAfterR_ = 1.0;
220              btranAverageAfterL_ = 1.0;
221            }
222          }
223          // scale back
224         
225          ftranCountInput_ *= 0.8;
226          ftranCountAfterL_ *= 0.8;
227          ftranCountAfterR_ *= 0.8;
228          ftranCountAfterU_ *= 0.8;
229          btranCountInput_ *= 0.8;
230          btranCountAfterU_ *= 0.8;
231          btranCountAfterR_ *= 0.8;
232          btranCountAfterL_ *= 0.8;
233        }
234      } else if (status_ == -1&&(solveType==0||solveType==2)) {
235        // This needs redoing as it was merged coding - does not need array
236        int numberTotal = numberRows+numberColumns;
237        int * isBasic = new int [numberTotal];
238        int * rowIsBasic = isBasic+numberColumns;
239        int * columnIsBasic = isBasic;
240        for (i=0;i<numberTotal;i++) 
241          isBasic[i]=-1;
242        for (i=0;i<numberRowBasic;i++) {
243          int iRow = pivotVariable[i];
244          rowIsBasic[iRow]=1;
245        }
246        for (;i<numberBasic;i++) {
247          int iColumn = pivotVariable[i];
248          columnIsBasic[iColumn]=1;
249        }
250        numberBasic=0;
251        for (i=0;i<numberRows;i++) 
252          pivotVariable[i]=-1;
253        // mark as basic or non basic
254        for (i=0;i<numberRows;i++) {
255          if (rowIsBasic[i]>=0) {
256            if (pivotColumn_[numberBasic]>=0) 
257              rowIsBasic[i]=pivotColumn_[numberBasic];
258            else
259              rowIsBasic[i]=-1;
260            numberBasic++;
261          }
262        }
263        for (i=0;i<numberColumns;i++) {
264          if (columnIsBasic[i]>=0) {
265            if (pivotColumn_[numberBasic]>=0) 
266              columnIsBasic[i]=pivotColumn_[numberBasic];
267            else
268              columnIsBasic[i]=-1;
269            numberBasic++;
270          }
271        }
272        // leave pivotVariable in useful form for cleaning basis
273        int * pivotVariable = model->pivotVariable();
274        for (i=0;i<numberRows;i++) {
275          pivotVariable[i]=-1;
276        }
277        for (i=0;i<numberRows;i++) {
278          if (model->getRowStatus(i) == ClpSimplex::basic) {
279            int iPivot = rowIsBasic[i];
280            if (iPivot>=0) 
281              pivotVariable[iPivot]=i+numberColumns;
282          }
283        }
284        for (i=0;i<numberColumns;i++) {
285          if (model->getColumnStatus(i) == ClpSimplex::basic) {
286            int iPivot = columnIsBasic[i];
287            if (iPivot>=0) 
288              pivotVariable[iPivot]=i;
289          }
290        }
291        delete [] isBasic;
292        double * columnLower = model->lowerRegion();
293        double * columnUpper = model->upperRegion();
294        double * columnActivity = model->solutionRegion();
295        double * rowLower = model->lowerRegion(0);
296        double * rowUpper = model->upperRegion(0);
297        double * rowActivity = model->solutionRegion(0);
298        //redo basis - first take ALL columns out
299        int iColumn;
300        double largeValue = model->largeValue();
301        for (iColumn=0;iColumn<numberColumns;iColumn++) {
302          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
303            // take out
304            if (!valuesPass) {
305              double lower=columnLower[iColumn];
306              double upper=columnUpper[iColumn];
307              double value=columnActivity[iColumn];
308              if (lower>-largeValue||upper<largeValue) {
309                if (fabs(value-lower)<fabs(value-upper)) {
310                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
311                  columnActivity[iColumn]=lower;
312                } else {
313                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
314                  columnActivity[iColumn]=upper;
315                }
316              } else {
317                model->setColumnStatus(iColumn,ClpSimplex::isFree);
318              }
319            } else {
320              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
321            }
322          }
323        }
324        int iRow;
325        for (iRow=0;iRow<numberRows;iRow++) {
326          int iSequence=pivotVariable[iRow];
327          if (iSequence>=0) {
328            // basic
329            if (iSequence>=numberColumns) {
330              // slack in - leave
331              //assert (iSequence-numberColumns==iRow);
332            } else {
333              // put back structural
334              model->setColumnStatus(iSequence,ClpSimplex::basic);
335            }
336          } else {
337            // put in slack
338            model->setRowStatus(iRow,ClpSimplex::basic);
339          }
340        }
341        // signal repeat
342        status_=-99;
343        // set fixed if they are
344        for (iRow=0;iRow<numberRows;iRow++) {
345          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
346            if (rowLower[iRow]==rowUpper[iRow]) {
347              rowActivity[iRow]=rowLower[iRow];
348              model->setRowStatus(iRow,ClpSimplex::isFixed);
349            }
350          }
351        }
352        for (iColumn=0;iColumn<numberColumns;iColumn++) {
353          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
354            if (columnLower[iColumn]==columnUpper[iColumn]) {
355              columnActivity[iColumn]=columnLower[iColumn];
356              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
357            }
358          }
359        }
360      } 
361    }
362  } else {
363    // network - fake factorization - do nothing
364    status_=0;
365  }
366
367  return status_;
368}
369/* Replaces one Column to basis,
370   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
371   If checkBeforeModifying is true will do all accuracy checks
372   before modifying factorization.  Whether to set this depends on
373   speed considerations.  You could just do this on first iteration
374   after factorization and thereafter re-factorize
375   partial update already in U */
376int 
377ClpFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
378                      int pivotRow,
379                      double pivotCheck ,
380                      bool checkBeforeModifying)
381{
382  if (!networkBasis_) {
383    return CoinFactorization::replaceColumn(regionSparse,
384                                            pivotRow,
385                                            pivotCheck,
386                                            checkBeforeModifying);
387  } else {
388    if (doCheck) {
389      int returnCode = CoinFactorization::replaceColumn(regionSparse,
390                                                        pivotRow,
391                                                        pivotCheck,
392                                                        checkBeforeModifying);
393      networkBasis_->replaceColumn(regionSparse,
394                                   pivotRow);
395      return returnCode;
396    } else {
397      // increase number of pivots
398      numberPivots_++;
399      return networkBasis_->replaceColumn(regionSparse,
400                                   pivotRow);
401    }
402  }
403}
404
405/* Updates one column (FTRAN) from region2
406   number returned is negative if no room
407   region1 starts as zero and is zero at end */
408int 
409ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
410                                   CoinIndexedVector * regionSparse2)
411{
412#ifdef CLP_DEBUG
413  regionSparse->checkClear();
414#endif
415  if (!networkBasis_) {
416    collectStatistics_ = true;
417    int returnValue= CoinFactorization::updateColumnFT(regionSparse,
418                                                       regionSparse2);
419    collectStatistics_ = false;
420    return returnValue;
421  } else {
422#ifdef CHECK_NETWORK
423    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
424    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
425                                                       regionSparse2);
426    networkBasis_->updateColumn(regionSparse,save);
427    int i;
428    double * array = regionSparse2->denseVector();
429    double * array2 = save->denseVector();
430    for (i=0;i<numberRows_;i++) {
431      double value1 = array[i];
432      double value2 = array2[i];
433      assert (value1==value2);
434    }
435    delete save;
436    return returnCode;
437#else
438    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
439    return 1;
440#endif
441  }
442}
443/* Updates one column (FTRAN) from region2
444   number returned is negative if no room
445   region1 starts as zero and is zero at end */
446int 
447ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
448                                 CoinIndexedVector * regionSparse2,
449                                 bool noPermute) const
450{
451#ifdef CLP_DEBUG
452  regionSparse->checkClear();
453#endif
454  if (!networkBasis_) {
455    collectStatistics_ = true;
456    int returnValue= CoinFactorization::updateColumn(regionSparse,
457                                                     regionSparse2,
458                                                     noPermute);
459    collectStatistics_ = false;
460    return returnValue;
461  } else {
462#ifdef CHECK_NETWORK
463    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
464    int returnCode = CoinFactorization::updateColumn(regionSparse,
465                                                     regionSparse2,
466                                                     noPermute);
467    networkBasis_->updateColumn(regionSparse,save);
468    int i;
469    double * array = regionSparse2->denseVector();
470    double * array2 = save->denseVector();
471    for (i=0;i<numberRows_;i++) {
472      double value1 = array[i];
473      double value2 = array2[i];
474      assert (value1==value2);
475    }
476    delete save;
477    return returnCode;
478#else
479    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
480    return 1;
481#endif
482  }
483}
484/* Updates one column (BTRAN) from region2
485   region1 starts as zero and is zero at end */
486int 
487ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
488                          CoinIndexedVector * regionSparse2) const
489{
490  if (!networkBasis_) {
491    collectStatistics_ = true;
492    return CoinFactorization::updateColumnTranspose(regionSparse,
493                                                    regionSparse2);
494    collectStatistics_ = false;
495  } else {
496#ifdef CHECK_NETWORK
497      CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
498      int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
499                                                                regionSparse2);
500      networkBasis_->updateColumnTranspose(regionSparse,save);
501      int i;
502      double * array = regionSparse2->denseVector();
503      double * array2 = save->denseVector();
504      for (i=0;i<numberRows_;i++) {
505        double value1 = array[i];
506        double value2 = array2[i];
507        assert (value1==value2);
508      }
509      delete save;
510      return returnCode;
511#else
512      return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
513#endif
514  }
515}
516/* makes a row copy of L for speed and to allow very sparse problems */
517void 
518ClpFactorization::goSparse()
519{
520  if (!networkBasis_) 
521    CoinFactorization::goSparse();
522}
523// Cleans up i.e. gets rid of network basis
524void 
525ClpFactorization::cleanUp()
526{
527  delete networkBasis_;
528  networkBasis_=NULL;
529  resetStatistics();
530}
531/// Says whether to redo pivot order
532bool 
533ClpFactorization::needToReorder() const
534{
535#ifdef CHECK_NETWORK
536  return true;
537#endif
538  if (!networkBasis_)
539    return true;
540  else
541    return false;
542}
Note: See TracBrowser for help on using the repository browser.