source: trunk/ClpFactorization.cpp @ 225

Last change on this file since 225 was 225, checked in by forrest, 16 years ago

This should break everything

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