source: trunk/ClpFactorization.cpp @ 454

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

Fix Network code

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