source: trunk/ClpFactorization.cpp @ 284

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

Wssmp cholesky and stuff for gub

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