source: trunk/ClpFactorization.cpp @ 461

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

Trying to make faster

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