source: trunk/ClpFactorization.cpp @ 338

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

NO_RTTI coding

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.4 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  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      // Move pivot variables across if they look good
96      int * pivotTemp = model->rowArray(0)->getIndices();
97      assert (!model->rowArray(0)->getNumElements());
98      if (!matrix->rhsOffset(model)) {
99        // Seems to prefer things in order so quickest
100        // way is to go though like this
101        for (i=0;i<numberRows;i++) {
102          if (model->getRowStatus(i) == ClpSimplex::basic) 
103            pivotTemp[numberBasic++]=i;
104        }
105        numberRowBasic=numberBasic;
106        /* Put column basic variables into pivotVariable
107           This is done by ClpMatrixBase to allow override for gub
108        */
109        matrix->generalExpanded(model,0,numberBasic);
110      } else {
111        // Long matrix - do a different way
112        bool fullSearch=false;
113        for (i=0;i<numberRows;i++) {
114          int iPivot = pivotVariable[i];
115          if (iPivot>=numberColumns) {
116            pivotTemp[numberBasic++]=iPivot-numberColumns;
117          }
118        }
119        numberRowBasic=numberBasic;
120        for (i=0;i<numberRows;i++) {
121          int iPivot = pivotVariable[i];
122          if (iPivot<numberColumns) {
123            if (iPivot>=0) {
124              pivotTemp[numberBasic++]=iPivot;
125            } else {
126              // not full basis
127              fullSearch=true;
128              break;
129            }
130          }
131        }
132        if (fullSearch) {
133          // do slow way
134          numberBasic=0;
135          for (i=0;i<numberRows;i++) {
136            if (model->getRowStatus(i) == ClpSimplex::basic) 
137              pivotTemp[numberBasic++]=i;
138          }
139          numberRowBasic=numberBasic;
140          /* Put column basic variables into pivotVariable
141             This is done by ClpMatrixBase to allow override for gub
142          */
143          matrix->generalExpanded(model,0,numberBasic);
144        }
145      }
146      assert (numberBasic<=model->maximumBasic());
147      // see if matrix a network
148#ifndef NO_RTTI
149      ClpNetworkMatrix* networkMatrix =
150        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
151#else
152      ClpNetworkMatrix* networkMatrix = NULL;
153      if (model->clpMatrix()->type()==11)
154        networkMatrix = 
155        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
156#endif
157      // If network - still allow ordinary factorization first time for laziness
158      int saveMaximumPivots = maximumPivots();
159      delete networkBasis_;
160      networkBasis_ = NULL;
161      if (networkMatrix&&!doCheck)
162        maximumPivots(1);
163      while (status_==-99) {
164        // maybe for speed will be better to leave as many regions as possible
165        gutsOfDestructor();
166        gutsOfInitialize(2);
167        CoinBigIndex numberElements=numberRowBasic;
168
169        // compute how much in basis
170
171        int i;
172        // can change for gub
173        int numberColumnBasic = numberBasic-numberRowBasic;
174
175        numberElements +=matrix->fillBasis(model,
176                                           pivotTemp+numberRowBasic, 
177                                           numberRowBasic,
178                                           numberColumnBasic,
179                                           NULL,NULL,NULL);
180        // and recompute as network side say different
181        if (model->numberIterations())
182          numberRowBasic = numberBasic - numberColumnBasic;
183        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
184        // If iteration not zero then may be compressed
185        getAreas ( !model->numberIterations() ? numberRows : numberBasic, 
186                   numberRowBasic+numberColumnBasic, numberElements,
187                   2 * numberElements );
188        //fill
189        //copy
190        numberElements=numberRowBasic;
191        for (i=0;i<numberRowBasic;i++) {
192          int iRow = pivotTemp[i];
193          indexRowU_[i]=iRow;
194          indexColumnU_[i]=i;
195          elementU_[i]=slackValue_;
196        }
197        // can change for gub so redo
198        numberColumnBasic = numberBasic-numberRowBasic;
199        numberElements +=matrix->fillBasis(model, 
200                                           pivotTemp+numberRowBasic, 
201                                           numberRowBasic, 
202                                           numberColumnBasic,
203                                           indexRowU_+numberElements, 
204                                           indexColumnU_+numberElements,
205                                           elementU_+numberElements);
206#if 0
207        {
208          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
209          for (int i=0;i<numberElements;i++)
210            printf("row %d col %d value %g\n",indexRowU_[i],indexColumnU_[i],
211                   elementU_[i]);
212        }
213#endif
214        // recompute number basic
215        numberBasic = numberRowBasic+numberColumnBasic;
216        lengthU_ = numberElements;
217
218        preProcess ( 0 );
219        factor (  );
220        if (status_==-99) {
221          // get more memory
222          areaFactor(2.0*areaFactor());
223        }
224      }
225      // If we get here status is 0 or -1
226      if (status_ == 0) {
227        // We may need to tamper with order and redo - e.g. network with side
228        int useNumberRows = numberRows;
229        // **** we will also need to add test in dual steepest to do
230        // as we do for network
231        matrix->generalExpanded(model,12,useNumberRows);
232        int * permuteBack = permuteBack_;
233        int * back = pivotColumnBack_;
234        //int * pivotTemp = pivotColumn_;
235        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
236        // Redo pivot order
237        for (i=0;i<numberRowBasic;i++) {
238          int k = pivotTemp[i];
239          // so rowIsBasic[k] would be permuteBack[back[i]]
240          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
241        }
242        for (;i<useNumberRows;i++) {
243          int k = pivotTemp[i];
244          // so rowIsBasic[k] would be permuteBack[back[i]]
245          pivotVariable[permuteBack[back[i]]]=k;
246        }
247        // Set up permutation vector
248        // these arrays start off as copies of permute
249        // (and we could use permute_ instead of pivotColumn (not back though))
250        ClpDisjointCopyN ( permute_, useNumberRows , pivotColumn_  );
251        ClpDisjointCopyN ( permuteBack_, useNumberRows , pivotColumnBack_  );
252        if (networkMatrix) {
253          maximumPivots(saveMaximumPivots);
254          // create network factorization
255          if (doCheck)
256            delete networkBasis_; // temp
257          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
258                                              pivotRegion_,
259                                              permuteBack_,
260                                              startColumnU_,
261                                              numberInColumn_,
262                                              indexRowU_,
263                                              elementU_);
264          // kill off arrays in ordinary factorization
265          if (!doCheck) {
266            gutsOfDestructor();
267            status_=0;
268#if 0
269            // but put back permute arrays so odd things will work
270            int numberRows = model->numberRows();
271            pivotColumnBack_ = new int [numberRows];
272            permute_ = new int [numberRows];
273            int i;
274            for (i=0;i<numberRows;i++) {
275              pivotColumnBack_[i]=i;
276              permute_[i]=i;
277            }
278#endif
279          }
280        } else {
281          // See if worth going sparse and when
282          if (numberFtranCounts_>100) {
283            ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
284            ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
285            ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
286            assert (ftranCountInput_&&ftranCountAfterL_&&ftranCountAfterR_);
287            if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
288              btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
289              btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
290              btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
291            } else {
292              // we have not done any useful btrans (values pass?)
293              btranAverageAfterU_ = 1.0;
294              btranAverageAfterR_ = 1.0;
295              btranAverageAfterL_ = 1.0;
296            }
297          }
298          // scale back
299         
300          ftranCountInput_ *= 0.8;
301          ftranCountAfterL_ *= 0.8;
302          ftranCountAfterR_ *= 0.8;
303          ftranCountAfterU_ *= 0.8;
304          btranCountInput_ *= 0.8;
305          btranCountAfterU_ *= 0.8;
306          btranCountAfterR_ *= 0.8;
307          btranCountAfterL_ *= 0.8;
308        }
309      } else if (status_ == -1&&(solveType==0||solveType==2)) {
310        // This needs redoing as it was merged coding - does not need array
311        int numberTotal = numberRows+numberColumns;
312        int * isBasic = new int [numberTotal];
313        int * rowIsBasic = isBasic+numberColumns;
314        int * columnIsBasic = isBasic;
315        for (i=0;i<numberTotal;i++) 
316          isBasic[i]=-1;
317        for (i=0;i<numberRowBasic;i++) {
318          int iRow = pivotTemp[i];
319          rowIsBasic[iRow]=1;
320        }
321        for (;i<numberBasic;i++) {
322          int iColumn = pivotTemp[i];
323          columnIsBasic[iColumn]=1;
324        }
325        numberBasic=0;
326        for (i=0;i<numberRows;i++) 
327          pivotVariable[i]=-1;
328        // mark as basic or non basic
329        for (i=0;i<numberRows;i++) {
330          if (rowIsBasic[i]>=0) {
331            if (pivotColumn_[numberBasic]>=0) 
332              rowIsBasic[i]=pivotColumn_[numberBasic];
333            else
334              rowIsBasic[i]=-1;
335            numberBasic++;
336          }
337        }
338        for (i=0;i<numberColumns;i++) {
339          if (columnIsBasic[i]>=0) {
340            if (pivotColumn_[numberBasic]>=0) 
341              columnIsBasic[i]=pivotColumn_[numberBasic];
342            else
343              columnIsBasic[i]=-1;
344            numberBasic++;
345          }
346        }
347        // leave pivotVariable in useful form for cleaning basis
348        int * pivotVariable = model->pivotVariable();
349        for (i=0;i<numberRows;i++) {
350          pivotVariable[i]=-1;
351        }
352        for (i=0;i<numberRows;i++) {
353          if (model->getRowStatus(i) == ClpSimplex::basic) {
354            int iPivot = rowIsBasic[i];
355            if (iPivot>=0) 
356              pivotVariable[iPivot]=i+numberColumns;
357          }
358        }
359        for (i=0;i<numberColumns;i++) {
360          if (model->getColumnStatus(i) == ClpSimplex::basic) {
361            int iPivot = columnIsBasic[i];
362            if (iPivot>=0) 
363              pivotVariable[iPivot]=i;
364          }
365        }
366        delete [] isBasic;
367        double * columnLower = model->lowerRegion();
368        double * columnUpper = model->upperRegion();
369        double * columnActivity = model->solutionRegion();
370        double * rowLower = model->lowerRegion(0);
371        double * rowUpper = model->upperRegion(0);
372        double * rowActivity = model->solutionRegion(0);
373        //redo basis - first take ALL columns out
374        int iColumn;
375        double largeValue = model->largeValue();
376        for (iColumn=0;iColumn<numberColumns;iColumn++) {
377          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
378            // take out
379            if (!valuesPass) {
380              double lower=columnLower[iColumn];
381              double upper=columnUpper[iColumn];
382              double value=columnActivity[iColumn];
383              if (lower>-largeValue||upper<largeValue) {
384                if (fabs(value-lower)<fabs(value-upper)) {
385                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
386                  columnActivity[iColumn]=lower;
387                } else {
388                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
389                  columnActivity[iColumn]=upper;
390                }
391              } else {
392                model->setColumnStatus(iColumn,ClpSimplex::isFree);
393              }
394            } else {
395              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
396            }
397          }
398        }
399        int iRow;
400        for (iRow=0;iRow<numberRows;iRow++) {
401          int iSequence=pivotVariable[iRow];
402          if (iSequence>=0) {
403            // basic
404            if (iSequence>=numberColumns) {
405              // slack in - leave
406              //assert (iSequence-numberColumns==iRow);
407            } else {
408              // put back structural
409              model->setColumnStatus(iSequence,ClpSimplex::basic);
410            }
411          } else {
412            // put in slack
413            model->setRowStatus(iRow,ClpSimplex::basic);
414          }
415        }
416        // Put back any key variables for gub (status_ not touched)
417        matrix->generalExpanded(model,1,status_);
418        // signal repeat
419        status_=-99;
420        // set fixed if they are
421        for (iRow=0;iRow<numberRows;iRow++) {
422          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
423            if (rowLower[iRow]==rowUpper[iRow]) {
424              rowActivity[iRow]=rowLower[iRow];
425              model->setRowStatus(iRow,ClpSimplex::isFixed);
426            }
427          }
428        }
429        for (iColumn=0;iColumn<numberColumns;iColumn++) {
430          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
431            if (columnLower[iColumn]==columnUpper[iColumn]) {
432              columnActivity[iColumn]=columnLower[iColumn];
433              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
434            }
435          }
436        }
437      } 
438    }
439  } else {
440    // network - fake factorization - do nothing
441    status_=0;
442  }
443
444  if (!status_) {
445    // take out part if quadratic
446    if (model->algorithm()==2) {
447      ClpObjective * obj = model->objectiveAsObject();
448      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
449      assert (quadraticObj);
450      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
451      int numberXColumns = quadratic->getNumCols();
452      assert (numberXColumns<numberColumns);
453      int base = numberColumns-numberXColumns;
454      int * which = new int [numberXColumns];
455      int * pivotVariable = model->pivotVariable();
456      int * permute = pivotColumn();
457      int i;
458      int n=0;
459      for (i=0;i<numberRows;i++) {
460        int iSj = pivotVariable[i]-base;
461        if (iSj>=0&&iSj<numberXColumns) 
462          which[n++]=permute[i];
463      }
464      if (n)
465        emptyRows(n,which);
466      delete [] which;
467    }
468  }
469  return status_;
470}
471/* Replaces one Column to basis,
472   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
473   If checkBeforeModifying is true will do all accuracy checks
474   before modifying factorization.  Whether to set this depends on
475   speed considerations.  You could just do this on first iteration
476   after factorization and thereafter re-factorize
477   partial update already in U */
478int 
479ClpFactorization::replaceColumn ( const ClpSimplex * model, 
480                                  CoinIndexedVector * regionSparse,
481                                  CoinIndexedVector * tableauColumn,
482                                  int pivotRow,
483                                  double pivotCheck ,
484                                  bool checkBeforeModifying)
485{
486  if (!networkBasis_) {
487    // see if FT
488    if (doForrestTomlin_)
489      return CoinFactorization::replaceColumn(regionSparse,
490                                              pivotRow,
491                                              pivotCheck,
492                                              checkBeforeModifying);
493    else
494      return CoinFactorization::replaceColumnPFI(tableauColumn,
495                                              pivotRow,pivotCheck); // Note array
496
497  } else {
498    if (doCheck) {
499      int returnCode = CoinFactorization::replaceColumn(regionSparse,
500                                                        pivotRow,
501                                                        pivotCheck,
502                                                        checkBeforeModifying);
503      networkBasis_->replaceColumn(regionSparse,
504                                   pivotRow);
505      return returnCode;
506    } else {
507      // increase number of pivots
508      numberPivots_++;
509      return networkBasis_->replaceColumn(regionSparse,
510                                   pivotRow);
511    }
512  }
513}
514
515/* Updates one column (FTRAN) from region2
516   number returned is negative if no room
517   region1 starts as zero and is zero at end */
518int 
519ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
520                                   CoinIndexedVector * regionSparse2)
521{
522#ifdef CLP_DEBUG
523  regionSparse->checkClear();
524#endif
525  if (!networkBasis_) {
526    collectStatistics_ = true;
527    int returnValue= CoinFactorization::updateColumnFT(regionSparse,
528                                                       regionSparse2);
529    collectStatistics_ = false;
530    return returnValue;
531  } else {
532#ifdef CHECK_NETWORK
533    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
534    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
535                                                       regionSparse2);
536    networkBasis_->updateColumn(regionSparse,save);
537    int i;
538    double * array = regionSparse2->denseVector();
539    double * array2 = save->denseVector();
540    for (i=0;i<numberRows_;i++) {
541      double value1 = array[i];
542      double value2 = array2[i];
543      assert (value1==value2);
544    }
545    delete save;
546    return returnCode;
547#else
548    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
549    return 1;
550#endif
551  }
552}
553/* Updates one column (FTRAN) from region2
554   number returned is negative if no room
555   region1 starts as zero and is zero at end */
556int 
557ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
558                                 CoinIndexedVector * regionSparse2,
559                                 bool noPermute) const
560{
561#ifdef CLP_DEBUG
562  if (!noPermute)
563    regionSparse->checkClear();
564#endif
565  if (!networkBasis_) {
566    collectStatistics_ = true;
567    int returnValue= CoinFactorization::updateColumn(regionSparse,
568                                                     regionSparse2,
569                                                     noPermute);
570    collectStatistics_ = false;
571    return returnValue;
572  } else {
573#ifdef CHECK_NETWORK
574    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
575    int returnCode = CoinFactorization::updateColumn(regionSparse,
576                                                     regionSparse2,
577                                                     noPermute);
578    networkBasis_->updateColumn(regionSparse,save);
579    int i;
580    double * array = regionSparse2->denseVector();
581    double * array2 = save->denseVector();
582    for (i=0;i<numberRows_;i++) {
583      double value1 = array[i];
584      double value2 = array2[i];
585      assert (value1==value2);
586    }
587    delete save;
588    return returnCode;
589#else
590    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
591    return 1;
592#endif
593  }
594}
595/* Updates one column (BTRAN) from region2
596   region1 starts as zero and is zero at end */
597int 
598ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
599                          CoinIndexedVector * regionSparse2) const
600{
601  if (!networkBasis_) {
602    collectStatistics_ = true;
603    return CoinFactorization::updateColumnTranspose(regionSparse,
604                                                    regionSparse2);
605    collectStatistics_ = false;
606  } else {
607#ifdef CHECK_NETWORK
608      CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
609      int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
610                                                                regionSparse2);
611      networkBasis_->updateColumnTranspose(regionSparse,save);
612      int i;
613      double * array = regionSparse2->denseVector();
614      double * array2 = save->denseVector();
615      for (i=0;i<numberRows_;i++) {
616        double value1 = array[i];
617        double value2 = array2[i];
618        assert (value1==value2);
619      }
620      delete save;
621      return returnCode;
622#else
623      return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
624#endif
625  }
626}
627/* makes a row copy of L for speed and to allow very sparse problems */
628void 
629ClpFactorization::goSparse()
630{
631  if (!networkBasis_) 
632    CoinFactorization::goSparse();
633}
634// Cleans up i.e. gets rid of network basis
635void 
636ClpFactorization::cleanUp()
637{
638  delete networkBasis_;
639  networkBasis_=NULL;
640  resetStatistics();
641}
642/// Says whether to redo pivot order
643bool 
644ClpFactorization::needToReorder() const
645{
646#ifdef CHECK_NETWORK
647  return true;
648#endif
649  if (!networkBasis_)
650    return true;
651  else
652    return false;
653}
Note: See TracBrowser for help on using the repository browser.