source: trunk/Clp/src/ClpFactorization.cpp @ 754

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.5 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#ifndef SLIM_CLP
7#include "ClpQuadraticObjective.hpp"
8#endif
9#include "CoinHelperFunctions.hpp"
10#include "CoinIndexedVector.hpp"
11#include "ClpSimplex.hpp"
12#include "ClpMatrixBase.hpp"
13#ifndef SLIM_CLP
14#include "ClpNetworkBasis.hpp"
15#include "ClpNetworkMatrix.hpp"
16//#define CHECK_NETWORK
17#ifdef CHECK_NETWORK
18const static bool doCheck=true;
19#else
20const static bool doCheck=false;
21#endif
22#endif
23//#############################################################################
24// Constructors / Destructor / Assignment
25//#############################################################################
26
27//-------------------------------------------------------------------
28// Default Constructor
29//-------------------------------------------------------------------
30ClpFactorization::ClpFactorization () :
31   CoinFactorization() 
32{
33#ifndef SLIM_CLP
34  networkBasis_ = NULL;
35#endif
36}
37
38//-------------------------------------------------------------------
39// Copy constructor
40//-------------------------------------------------------------------
41ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
42   CoinFactorization(rhs) 
43{
44#ifndef SLIM_CLP
45  if (rhs.networkBasis_)
46    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
47  else
48    networkBasis_=NULL;
49#endif
50}
51
52ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
53   CoinFactorization(rhs) 
54{
55#ifndef SLIM_CLP
56  networkBasis_=NULL;
57#endif
58}
59
60//-------------------------------------------------------------------
61// Destructor
62//-------------------------------------------------------------------
63ClpFactorization::~ClpFactorization () 
64{
65#ifndef SLIM_CLP
66  delete networkBasis_;
67#endif
68}
69
70//----------------------------------------------------------------
71// Assignment operator
72//-------------------------------------------------------------------
73ClpFactorization &
74ClpFactorization::operator=(const ClpFactorization& rhs)
75{
76  if (this != &rhs) {
77    CoinFactorization::operator=(rhs);
78#ifndef SLIM_CLP
79    delete networkBasis_;
80    if (rhs.networkBasis_)
81      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
82    else
83      networkBasis_=NULL;
84#endif
85  }
86  return *this;
87}
88#if 0
89static unsigned int saveList[10000];
90int numberSave=-1;
91inline bool isDense(int i) {
92  return ((saveList[i>>5]>>(i&31))&1)!=0;
93}
94inline void setDense(int i) {
95  unsigned int & value = saveList[i>>5];
96  int bit = i&31;
97  value |= (1<<bit);
98}
99#endif
100int 
101ClpFactorization::factorize ( ClpSimplex * model,
102                              int solveType, bool valuesPass)
103{
104  ClpMatrixBase * matrix = model->clpMatrix(); 
105  int numberRows = model->numberRows();
106  int numberColumns = model->numberColumns();
107  if (!numberRows)
108    return 0;
109  // If too many compressions increase area
110  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
111    areaFactor_ *= 1.1;
112  }
113  //int numberPivots=numberPivots_;
114#if 0
115  if (model->algorithm()>0)
116    numberSave=-1;
117#endif
118#ifndef SLIM_CLP
119  if (!networkBasis_||doCheck) {
120#endif
121    status_=-99;
122    int * pivotVariable = model->pivotVariable();
123    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
124    while (status_<-98) {
125     
126      int i;
127      int numberBasic=0;
128      int numberRowBasic;
129      // Move pivot variables across if they look good
130      int * pivotTemp = model->rowArray(0)->getIndices();
131      assert (!model->rowArray(0)->getNumElements());
132      if (!matrix->rhsOffset(model)) {
133#if 0
134        if (numberSave>0) {
135          int nStill=0;
136          int nAtBound=0;
137          int nZeroDual=0;
138          CoinIndexedVector * array = model->rowArray(3);
139          CoinIndexedVector * objArray = model->columnArray(1);
140          array->clear();
141          objArray->clear();
142          double * cost = model->costRegion();
143          double tolerance = model->primalTolerance();
144          double offset=0.0;
145          for (i=0;i<numberRows;i++) {
146            int iPivot = pivotVariable[i];
147            if (iPivot<numberColumns&&isDense(iPivot)) {
148              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
149                nStill++;
150                double value=model->solutionRegion()[iPivot];
151                double dual = model->dualRowSolution()[i];
152                double lower=model->lowerRegion()[iPivot];
153                double upper=model->upperRegion()[iPivot];
154                ClpSimplex::Status status;
155                if (fabs(value-lower)<tolerance) {
156                  status=ClpSimplex::atLowerBound;
157                  nAtBound++;
158                } else if (fabs(value-upper)<tolerance) {
159                  nAtBound++;
160                  status=ClpSimplex::atUpperBound;
161                } else if (value>lower&&value<upper) {
162                  status=ClpSimplex::superBasic;
163                } else {
164                  status=ClpSimplex::basic;
165                }
166                if (status!=ClpSimplex::basic) {
167                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
168                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
169                    model->setRowStatus(i,ClpSimplex::basic);
170                    pivotVariable[i]=i+numberColumns;
171                    model->dualRowSolution()[i]=0.0;
172                    model->djRegion(0)[i]=0.0;
173                    array->add(i,dual);
174                    offset += dual*model->solutionRegion(0)[i];
175                  }
176                }
177                if (fabs(dual)<1.0e-5)
178                  nZeroDual++;
179              }
180            }
181          }
182          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
183                 numberSave,nStill,nAtBound,nZeroDual,offset);
184          if (array->getNumElements()) {
185            // modify costs
186            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
187                                               objArray);
188            array->clear();
189            int n=objArray->getNumElements();
190            int * indices = objArray->getIndices();
191            double * elements = objArray->denseVector();
192            for (i=0;i<n;i++) {
193              int iColumn = indices[i];
194              cost[iColumn] -= elements[iColumn];
195              elements[iColumn]=0.0;
196            }
197            objArray->setNumElements(0);
198          }
199        }
200#endif
201        // Seems to prefer things in order so quickest
202        // way is to go though like this
203        for (i=0;i<numberRows;i++) {
204          if (model->getRowStatus(i) == ClpSimplex::basic) 
205            pivotTemp[numberBasic++]=i;
206        }
207        numberRowBasic=numberBasic;
208        /* Put column basic variables into pivotVariable
209           This is done by ClpMatrixBase to allow override for gub
210        */
211        matrix->generalExpanded(model,0,numberBasic);
212      } else {
213        // Long matrix - do a different way
214        bool fullSearch=false;
215        for (i=0;i<numberRows;i++) {
216          int iPivot = pivotVariable[i];
217          if (iPivot>=numberColumns) {
218            pivotTemp[numberBasic++]=iPivot-numberColumns;
219          }
220        }
221        numberRowBasic=numberBasic;
222        for (i=0;i<numberRows;i++) {
223          int iPivot = pivotVariable[i];
224          if (iPivot<numberColumns) {
225            if (iPivot>=0) {
226              pivotTemp[numberBasic++]=iPivot;
227            } else {
228              // not full basis
229              fullSearch=true;
230              break;
231            }
232          }
233        }
234        if (fullSearch) {
235          // do slow way
236          numberBasic=0;
237          for (i=0;i<numberRows;i++) {
238            if (model->getRowStatus(i) == ClpSimplex::basic) 
239              pivotTemp[numberBasic++]=i;
240          }
241          numberRowBasic=numberBasic;
242          /* Put column basic variables into pivotVariable
243             This is done by ClpMatrixBase to allow override for gub
244          */
245          matrix->generalExpanded(model,0,numberBasic);
246        }
247      }
248      if (numberBasic>model->maximumBasic()) {
249#ifndef NDEBUG
250        printf("%d basic - should only be %d\n",
251               numberBasic,numberRows);
252#endif
253        // Take out some
254        numberBasic=numberRowBasic;
255        for (int i=0;i<numberColumns;i++) {
256          if (model->getColumnStatus(i) == ClpSimplex::basic) {
257            if (numberBasic<numberRows)
258              numberBasic++;
259            else
260              model->setColumnStatus(i,ClpSimplex::superBasic);
261          }
262        }
263        numberBasic=numberRowBasic;
264        matrix->generalExpanded(model,0,numberBasic);
265      }
266#ifndef SLIM_CLP
267      // see if matrix a network
268#ifndef NO_RTTI
269      ClpNetworkMatrix* networkMatrix =
270        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
271#else
272      ClpNetworkMatrix* networkMatrix = NULL;
273      if (model->clpMatrix()->type()==11)
274        networkMatrix = 
275        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
276#endif
277      // If network - still allow ordinary factorization first time for laziness
278      if (networkMatrix)
279        biasLU_=0; // All to U if network
280      int saveMaximumPivots = maximumPivots();
281      delete networkBasis_;
282      networkBasis_ = NULL;
283      if (networkMatrix&&!doCheck)
284        maximumPivots(1);
285#endif
286      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
287      while (status_==-99) {
288        // maybe for speed will be better to leave as many regions as possible
289        gutsOfDestructor();
290        gutsOfInitialize(2);
291        CoinBigIndex numberElements=numberRowBasic;
292
293        // compute how much in basis
294
295        int i;
296        // can change for gub
297        int numberColumnBasic = numberBasic-numberRowBasic;
298
299        numberElements +=matrix->countBasis(model,
300                                           pivotTemp+numberRowBasic, 
301                                           numberRowBasic,
302                                            numberColumnBasic);
303        // and recompute as network side say different
304        if (model->numberIterations())
305          numberRowBasic = numberBasic - numberColumnBasic;
306        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
307#if 0
308        // If iteration not zero then may be compressed
309        getAreas ( !model->numberIterations() ? numberRows : numberBasic,
310                   numberRowBasic+numberColumnBasic, numberElements,
311                   2 * numberElements );
312#else
313        getAreas ( numberRows,
314                   numberRowBasic+numberColumnBasic, numberElements,
315                   2 * numberElements );
316#endif
317        //fill
318        // Fill in counts so we can skip part of preProcess
319        CoinZeroN ( numberInRow_, numberRows_ + 1 );
320        CoinZeroN ( numberInColumn_, maximumColumnsExtra_ + 1 );
321        for (i=0;i<numberRowBasic;i++) {
322          int iRow = pivotTemp[i];
323          indexRowU_[i]=iRow;
324          startColumnU_[i]=i;
325          elementU_[i]=slackValue_;
326          numberInRow_[iRow]=1;
327          numberInColumn_[i]=1;
328        }
329        startColumnU_[numberRowBasic]=numberRowBasic;
330        // can change for gub so redo
331        numberColumnBasic = numberBasic-numberRowBasic;
332        matrix->fillBasis(model, 
333                          pivotTemp+numberRowBasic, 
334                          numberColumnBasic,
335                          indexRowU_, 
336                          startColumnU_+numberRowBasic,
337                          numberInRow_,
338                          numberInColumn_+numberRowBasic,
339                          elementU_);
340#if 0
341        {
342          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
343          for (int i=0;i<numberElements;i++)
344            printf("row %d col %d value %g\n",indexRowU_[i],indexColumnU_[i],
345                   elementU_[i]);
346        }
347#endif
348        // recompute number basic
349        numberBasic = numberRowBasic+numberColumnBasic;
350        if (numberBasic) 
351          numberElements = startColumnU_[numberBasic-1]
352            +numberInColumn_[numberBasic-1];
353        else
354          numberElements=0;
355        lengthU_ = numberElements;
356        //saveFactorization("dump.d");
357        if (biasLU_>=3||numberRows_!=numberColumns_)
358          preProcess ( 2 );
359        else
360          preProcess ( 3 ); // no row copy
361        factor (  );
362        if (status_==-99) {
363          // get more memory
364          areaFactor(2.0*areaFactor());
365        } else if (status_==-1&&model->numberIterations()==0&&
366                   denseThreshold_) {
367          // Round again without dense
368          denseThreshold_=0;
369          status_=-99;
370        }
371      }
372      // If we get here status is 0 or -1
373      if (status_ == 0) {
374        // We may need to tamper with order and redo - e.g. network with side
375        int useNumberRows = numberRows;
376        // **** we will also need to add test in dual steepest to do
377        // as we do for network
378        matrix->generalExpanded(model,12,useNumberRows);
379        int * permuteBack = permuteBack_;
380        int * back = pivotColumnBack_;
381        //int * pivotTemp = pivotColumn_;
382        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
383        // Redo pivot order
384        for (i=0;i<numberRowBasic;i++) {
385          int k = pivotTemp[i];
386          // so rowIsBasic[k] would be permuteBack[back[i]]
387          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
388        }
389        for (;i<useNumberRows;i++) {
390          int k = pivotTemp[i];
391          // so rowIsBasic[k] would be permuteBack[back[i]]
392          pivotVariable[permuteBack[back[i]]]=k;
393        }
394#if 0
395        if (numberSave>=0) {
396          numberSave=numberDense_;
397          memset(saveList,0,((numberRows_+31)>>5)*sizeof(int));
398          for (i=numberRows_-numberSave;i<numberRows_;i++) {
399            int k=pivotTemp[pivotColumn_[i]];
400            setDense(k);
401          }
402        }
403#endif
404        // Set up permutation vector
405        // these arrays start off as copies of permute
406        // (and we could use permute_ instead of pivotColumn (not back though))
407        ClpDisjointCopyN ( permute_, useNumberRows , pivotColumn_  );
408        ClpDisjointCopyN ( permuteBack_, useNumberRows , pivotColumnBack_  );
409#ifndef SLIM_CLP
410        if (networkMatrix) {
411          maximumPivots(saveMaximumPivots);
412          // create network factorization
413          if (doCheck)
414            delete networkBasis_; // temp
415          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
416                                              pivotRegion_,
417                                              permuteBack_,
418                                              startColumnU_,
419                                              numberInColumn_,
420                                              indexRowU_,
421                                              elementU_);
422          // kill off arrays in ordinary factorization
423          if (!doCheck) {
424            gutsOfDestructor();
425            // but make sure numberRows_ set
426            numberRows_ = model->numberRows();
427            status_=0;
428#if 0
429            // but put back permute arrays so odd things will work
430            int numberRows = model->numberRows();
431            pivotColumnBack_ = new int [numberRows];
432            permute_ = new int [numberRows];
433            int i;
434            for (i=0;i<numberRows;i++) {
435              pivotColumnBack_[i]=i;
436              permute_[i]=i;
437            }
438#endif
439          }
440        } else {
441#endif
442          // See if worth going sparse and when
443          if (numberFtranCounts_>100) {
444            ftranCountInput_= CoinMax(ftranCountInput_,1.0);
445            ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,1.0);
446            ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,1.0);
447            ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,1.0);
448            if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
449              btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,1.0);
450              btranAverageAfterR_ = CoinMax(btranCountAfterR_/btranCountAfterU_,1.0);
451              btranAverageAfterL_ = CoinMax(btranCountAfterL_/btranCountAfterR_,1.0);
452            } else {
453              // we have not done any useful btrans (values pass?)
454              btranAverageAfterU_ = 1.0;
455              btranAverageAfterR_ = 1.0;
456              btranAverageAfterL_ = 1.0;
457            }
458          }
459          // scale back
460         
461          ftranCountInput_ *= 0.8;
462          ftranCountAfterL_ *= 0.8;
463          ftranCountAfterR_ *= 0.8;
464          ftranCountAfterU_ *= 0.8;
465          btranCountInput_ *= 0.8;
466          btranCountAfterU_ *= 0.8;
467          btranCountAfterR_ *= 0.8;
468          btranCountAfterL_ *= 0.8;
469#ifndef SLIM_CLP
470        }
471#endif
472      } else if (status_ == -1&&(solveType==0||solveType==2)) {
473        // This needs redoing as it was merged coding - does not need array
474        int numberTotal = numberRows+numberColumns;
475        int * isBasic = new int [numberTotal];
476        int * rowIsBasic = isBasic+numberColumns;
477        int * columnIsBasic = isBasic;
478        for (i=0;i<numberTotal;i++) 
479          isBasic[i]=-1;
480        for (i=0;i<numberRowBasic;i++) {
481          int iRow = pivotTemp[i];
482          rowIsBasic[iRow]=1;
483        }
484        for (;i<numberBasic;i++) {
485          int iColumn = pivotTemp[i];
486          columnIsBasic[iColumn]=1;
487        }
488        numberBasic=0;
489        for (i=0;i<numberRows;i++) 
490          pivotVariable[i]=-1;
491        // mark as basic or non basic
492        for (i=0;i<numberRows;i++) {
493          if (rowIsBasic[i]>=0) {
494            if (pivotColumn_[numberBasic]>=0) {
495              rowIsBasic[i]=pivotColumn_[numberBasic];
496            } else {
497              rowIsBasic[i]=-1;
498              model->setRowStatus(i,ClpSimplex::superBasic);
499            }
500            numberBasic++;
501          }
502        }
503        for (i=0;i<numberColumns;i++) {
504          if (columnIsBasic[i]>=0) {
505            if (pivotColumn_[numberBasic]>=0) 
506              columnIsBasic[i]=pivotColumn_[numberBasic];
507            else
508              columnIsBasic[i]=-1;
509            numberBasic++;
510          }
511        }
512        // leave pivotVariable in useful form for cleaning basis
513        int * pivotVariable = model->pivotVariable();
514        for (i=0;i<numberRows;i++) {
515          pivotVariable[i]=-1;
516        }
517
518        for (i=0;i<numberRows;i++) {
519          if (model->getRowStatus(i) == ClpSimplex::basic) {
520            int iPivot = rowIsBasic[i];
521            if (iPivot>=0) 
522              pivotVariable[iPivot]=i+numberColumns;
523          }
524        }
525        for (i=0;i<numberColumns;i++) {
526          if (model->getColumnStatus(i) == ClpSimplex::basic) {
527            int iPivot = columnIsBasic[i];
528            if (iPivot>=0) 
529              pivotVariable[iPivot]=i;
530          }
531        }
532        delete [] isBasic;
533        double * columnLower = model->lowerRegion();
534        double * columnUpper = model->upperRegion();
535        double * columnActivity = model->solutionRegion();
536        double * rowLower = model->lowerRegion(0);
537        double * rowUpper = model->upperRegion(0);
538        double * rowActivity = model->solutionRegion(0);
539        //redo basis - first take ALL columns out
540        int iColumn;
541        double largeValue = model->largeValue();
542        for (iColumn=0;iColumn<numberColumns;iColumn++) {
543          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
544            // take out
545            if (!valuesPass) {
546              double lower=columnLower[iColumn];
547              double upper=columnUpper[iColumn];
548              double value=columnActivity[iColumn];
549              if (lower>-largeValue||upper<largeValue) {
550                if (fabs(value-lower)<fabs(value-upper)) {
551                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
552                  columnActivity[iColumn]=lower;
553                } else {
554                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
555                  columnActivity[iColumn]=upper;
556                }
557              } else {
558                model->setColumnStatus(iColumn,ClpSimplex::isFree);
559              }
560            } else {
561              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
562            }
563          }
564        }
565        int iRow;
566        for (iRow=0;iRow<numberRows;iRow++) {
567          int iSequence=pivotVariable[iRow];
568          if (iSequence>=0) {
569            // basic
570            if (iSequence>=numberColumns) {
571              // slack in - leave
572              //assert (iSequence-numberColumns==iRow);
573            } else {
574              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
575              // put back structural
576              model->setColumnStatus(iSequence,ClpSimplex::basic);
577            }
578          } else {
579            // put in slack
580            model->setRowStatus(iRow,ClpSimplex::basic);
581          }
582        }
583        // Put back any key variables for gub (status_ not touched)
584        matrix->generalExpanded(model,1,status_);
585        // signal repeat
586        status_=-99;
587        // set fixed if they are
588        for (iRow=0;iRow<numberRows;iRow++) {
589          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
590            if (rowLower[iRow]==rowUpper[iRow]) {
591              rowActivity[iRow]=rowLower[iRow];
592              model->setRowStatus(iRow,ClpSimplex::isFixed);
593            }
594          }
595        }
596        for (iColumn=0;iColumn<numberColumns;iColumn++) {
597          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
598            if (columnLower[iColumn]==columnUpper[iColumn]) {
599              columnActivity[iColumn]=columnLower[iColumn];
600              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
601            }
602          }
603        }
604      } 
605    }
606#ifndef SLIM_CLP
607  } else {
608    // network - fake factorization - do nothing
609    status_=0;
610  }
611#endif
612#ifndef SLIM_CLP
613  if (!status_) {
614    // take out part if quadratic
615    if (model->algorithm()==2) {
616      ClpObjective * obj = model->objectiveAsObject();
617      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
618      assert (quadraticObj);
619      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
620      int numberXColumns = quadratic->getNumCols();
621      assert (numberXColumns<numberColumns);
622      int base = numberColumns-numberXColumns;
623      int * which = new int [numberXColumns];
624      int * pivotVariable = model->pivotVariable();
625      int * permute = pivotColumn();
626      int i;
627      int n=0;
628      for (i=0;i<numberRows;i++) {
629        int iSj = pivotVariable[i]-base;
630        if (iSj>=0&&iSj<numberXColumns) 
631          which[n++]=permute[i];
632      }
633      if (n)
634        emptyRows(n,which);
635      delete [] which;
636    }
637  }
638#endif
639  return status_;
640}
641/* Replaces one Column to basis,
642   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
643   If checkBeforeModifying is true will do all accuracy checks
644   before modifying factorization.  Whether to set this depends on
645   speed considerations.  You could just do this on first iteration
646   after factorization and thereafter re-factorize
647   partial update already in U */
648int 
649ClpFactorization::replaceColumn ( const ClpSimplex * model, 
650                                  CoinIndexedVector * regionSparse,
651                                  CoinIndexedVector * tableauColumn,
652                                  int pivotRow,
653                                  double pivotCheck ,
654                                  bool checkBeforeModifying)
655{
656#ifndef SLIM_CLP
657  if (!networkBasis_) {
658#endif
659    // see if FT
660    if (doForrestTomlin_) {
661      int returnCode= CoinFactorization::replaceColumn(regionSparse,
662                                              pivotRow,
663                                              pivotCheck,
664                                              checkBeforeModifying);
665      return returnCode;
666    } else {
667      return CoinFactorization::replaceColumnPFI(tableauColumn,
668                                              pivotRow,pivotCheck); // Note array
669    }
670
671#ifndef SLIM_CLP
672  } else {
673    if (doCheck) {
674      int returnCode = CoinFactorization::replaceColumn(regionSparse,
675                                                        pivotRow,
676                                                        pivotCheck,
677                                                        checkBeforeModifying);
678      networkBasis_->replaceColumn(regionSparse,
679                                   pivotRow);
680      return returnCode;
681    } else {
682      // increase number of pivots
683      numberPivots_++;
684      return networkBasis_->replaceColumn(regionSparse,
685                                   pivotRow);
686    }
687  }
688#endif
689}
690
691/* Updates one column (FTRAN) from region2
692   number returned is negative if no room
693   region1 starts as zero and is zero at end */
694int 
695ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
696                                   CoinIndexedVector * regionSparse2)
697{
698#ifdef CLP_DEBUG
699  regionSparse->checkClear();
700#endif
701  if (!numberRows_)
702    return 0;
703#ifndef SLIM_CLP
704  if (!networkBasis_) {
705#endif
706    collectStatistics_ = true;
707    int returnValue= CoinFactorization::updateColumnFT(regionSparse,
708                                                       regionSparse2);
709    collectStatistics_ = false;
710    return returnValue;
711#ifndef SLIM_CLP
712  } else {
713#ifdef CHECK_NETWORK
714    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
715    double * check = new double[numberRows_];
716    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
717                                                       regionSparse2);
718    networkBasis_->updateColumn(regionSparse,save,-1);
719    int i;
720    double * array = regionSparse2->denseVector();
721    int * indices = regionSparse2->getIndices();
722    int n=regionSparse2->getNumElements();
723    memset(check,0,numberRows_*sizeof(double));
724    double * array2 = save->denseVector();
725    int * indices2 = save->getIndices();
726    int n2=save->getNumElements();
727    assert (n==n2);
728    if (save->packedMode()) {
729      for (i=0;i<n;i++) {
730        check[indices[i]]=array[i];
731      }
732      for (i=0;i<n;i++) {
733        double value2 = array2[i];
734        assert (check[indices2[i]]==value2);
735      }
736    } else {
737      for (i=0;i<numberRows_;i++) {
738        double value1 = array[i];
739        double value2 = array2[i];
740        assert (value1==value2);
741      }
742    }
743    delete save;
744    delete [] check;
745    return returnCode;
746#else
747    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
748    return 1;
749#endif
750  }
751#endif
752}
753/* Updates one column (FTRAN) from region2
754   number returned is negative if no room
755   region1 starts as zero and is zero at end */
756int 
757ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
758                                 CoinIndexedVector * regionSparse2,
759                                 bool noPermute) const
760{
761#ifdef CLP_DEBUG
762  if (!noPermute)
763    regionSparse->checkClear();
764#endif
765  if (!numberRows_)
766    return 0;
767#ifndef SLIM_CLP
768  if (!networkBasis_) {
769#endif
770    collectStatistics_ = true;
771    int returnValue= CoinFactorization::updateColumn(regionSparse,
772                                                     regionSparse2,
773                                                     noPermute);
774    collectStatistics_ = false;
775    return returnValue;
776#ifndef SLIM_CLP
777  } else {
778#ifdef CHECK_NETWORK
779    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
780    double * check = new double[numberRows_];
781    int returnCode = CoinFactorization::updateColumn(regionSparse,
782                                                     regionSparse2,
783                                                     noPermute);
784    networkBasis_->updateColumn(regionSparse,save,-1);
785    int i;
786    double * array = regionSparse2->denseVector();
787    int * indices = regionSparse2->getIndices();
788    int n=regionSparse2->getNumElements();
789    memset(check,0,numberRows_*sizeof(double));
790    double * array2 = save->denseVector();
791    int * indices2 = save->getIndices();
792    int n2=save->getNumElements();
793    assert (n==n2);
794    if (save->packedMode()) {
795      for (i=0;i<n;i++) {
796        check[indices[i]]=array[i];
797      }
798      for (i=0;i<n;i++) {
799        double value2 = array2[i];
800        assert (check[indices2[i]]==value2);
801      }
802    } else {
803      for (i=0;i<numberRows_;i++) {
804        double value1 = array[i];
805        double value2 = array2[i];
806        assert (value1==value2);
807      }
808    }
809    delete save;
810    delete [] check;
811    return returnCode;
812#else
813    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
814    return 1;
815#endif
816  }
817#endif
818}
819/* Updates one column (FTRAN) from region2
820   number returned is negative if no room
821   region1 starts as zero and is zero at end */
822int 
823ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
824                                 CoinIndexedVector * regionSparse2,
825                                 bool noPermute) const
826{
827  if (!noPermute)
828    regionSparse->checkClear();
829  if (!numberRows_)
830    return 0;
831  collectStatistics_ = false;
832  int returnValue= CoinFactorization::updateColumn(regionSparse,
833                                                   regionSparse2,
834                                                   noPermute);
835  return returnValue;
836}
837/* Updates one column (BTRAN) from region2
838   region1 starts as zero and is zero at end */
839int 
840ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
841                          CoinIndexedVector * regionSparse2) const
842{
843  if (!numberRows_)
844    return 0;
845#ifndef SLIM_CLP
846  if (!networkBasis_) {
847#endif
848    collectStatistics_ = true;
849    return CoinFactorization::updateColumnTranspose(regionSparse,
850                                                    regionSparse2);
851    collectStatistics_ = false;
852#ifndef SLIM_CLP
853  } else {
854#ifdef CHECK_NETWORK
855    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
856    double * check = new double[numberRows_];
857    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
858                                                              regionSparse2);
859    networkBasis_->updateColumnTranspose(regionSparse,save);
860    int i;
861    double * array = regionSparse2->denseVector();
862    int * indices = regionSparse2->getIndices();
863    int n=regionSparse2->getNumElements();
864    memset(check,0,numberRows_*sizeof(double));
865    double * array2 = save->denseVector();
866    int * indices2 = save->getIndices();
867    int n2=save->getNumElements();
868    assert (n==n2);
869    if (save->packedMode()) {
870      for (i=0;i<n;i++) {
871        check[indices[i]]=array[i];
872      }
873      for (i=0;i<n;i++) {
874        double value2 = array2[i];
875        assert (check[indices2[i]]==value2);
876      }
877    } else {
878      for (i=0;i<numberRows_;i++) {
879        double value1 = array[i];
880        double value2 = array2[i];
881        assert (value1==value2);
882      }
883    }
884    delete save;
885    delete [] check;
886    return returnCode;
887#else
888    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
889#endif
890  }
891#endif
892}
893/* makes a row copy of L for speed and to allow very sparse problems */
894void 
895ClpFactorization::goSparse()
896{
897#ifndef SLIM_CLP
898  if (!networkBasis_) 
899#endif
900    CoinFactorization::goSparse();
901}
902// Cleans up i.e. gets rid of network basis
903void 
904ClpFactorization::cleanUp()
905{
906#ifndef SLIM_CLP
907  delete networkBasis_;
908  networkBasis_=NULL;
909#endif
910  resetStatistics();
911}
912/// Says whether to redo pivot order
913bool 
914ClpFactorization::needToReorder() const
915{
916#ifdef CHECK_NETWORK
917  return true;
918#endif
919#ifndef SLIM_CLP
920  if (!networkBasis_)
921#endif
922    return true;
923#ifndef SLIM_CLP
924  else
925    return false;
926#endif
927}
928// Get weighted row list
929void
930ClpFactorization::getWeights(int * weights) const
931{
932#ifndef SLIM_CLP
933  if (networkBasis_) {
934    // Network - just unit
935    for (int i=0;i<numberRows_;i++) 
936      weights[i]=1;
937    return;
938  }
939#endif
940  int * permuteBack = pivotColumnBack_;
941  if (!startRowL_||!numberInRow_) {
942    int * temp = new int[numberRows_];
943    memset(temp,0,numberRows_*sizeof(int));
944    int i;
945    for (i=0;i<numberRows_;i++) {
946      // one for pivot
947      temp[i]++;
948      CoinBigIndex j;
949      for (j=startColumnU_[i];j<startColumnU_[i]+numberInColumn_[i];j++) {
950        int iRow=indexRowU_[j];
951        temp[iRow]++;
952      }
953    }
954    for (i=baseL_;i<baseL_+numberL_;i++) {
955      CoinBigIndex j;
956      for (j=startColumnL_[i];j<startColumnL_[i+1];j++) {
957        int iRow = indexRowL_[j];
958        temp[iRow]++;
959      }
960    }
961    for (i=0;i<numberRows_;i++) {
962      int number = temp[i];
963      int iPermute = permuteBack[i];
964      weights[iPermute]=number;
965    }
966    delete [] temp;
967  } else {
968    int i;
969    for (i=0;i<numberRows_;i++) {
970      int number = startRowL_[i+1]-startRowL_[i]+numberInRow_[i]+1;
971      //number = startRowL_[i+1]-startRowL_[i]+1;
972      //number = numberInRow_[i]+1;
973      int iPermute = permuteBack[i];
974      weights[iPermute]=number;
975    }
976  }
977}
Note: See TracBrowser for help on using the repository browser.