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

Last change on this file since 1371 was 1371, checked in by forrest, 11 years ago

changes to try and make faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 84.2 KB
Line 
1/* $Id: ClpFactorization.cpp 1371 2009-06-12 16:29:04Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6#include "ClpFactorization.hpp"
7#ifndef SLIM_CLP
8#include "ClpQuadraticObjective.hpp"
9#endif
10#include "CoinHelperFunctions.hpp"
11#include "CoinIndexedVector.hpp"
12#include "ClpSimplex.hpp"
13#include "ClpSimplexDual.hpp"
14#include "ClpMatrixBase.hpp"
15#ifndef SLIM_CLP
16#include "ClpNetworkBasis.hpp"
17#include "ClpNetworkMatrix.hpp"
18//#define CHECK_NETWORK
19#ifdef CHECK_NETWORK
20const static bool doCheck=true;
21#else
22const static bool doCheck=false;
23#endif
24#endif
25//#define CLP_FACTORIZATION_INSTRUMENT
26#ifdef CLP_FACTORIZATION_INSTRUMENT
27#include "CoinTime.hpp"
28double factorization_instrument(int type)
29{
30  static int times[10]={0,0,0,0,0,0,0,0,0,0};
31  static double startTime=0.0;
32  static double totalTimes [10]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
33  if (type<0) {
34    assert (!startTime);
35    startTime=CoinCpuTime();
36    return 0.0;
37  } else if (type>0) {
38    times[type]++;
39    double difference = CoinCpuTime()-startTime;
40    totalTimes[type] += difference;
41    startTime=0.0;
42    return difference;
43  } else {
44    // report
45    const char * types[10]=
46      {"","fac=rhs_etc","factorize","replace","update_FT",
47       "update","update_transpose","gosparse","getWeights!","update2_FT"};
48    double total=0.0;
49    for (int i=1;i<10;i++) {
50      if (times[i]) {
51        printf("%s was called %d times taking %g seconds\n",
52               types[i],times[i],totalTimes[i]);
53        total += totalTimes[i];
54        times[i]=0;
55        totalTimes[i]=0.0;
56      }
57    }
58    return total;
59  }
60}
61#endif
62//#############################################################################
63// Constructors / Destructor / Assignment
64//#############################################################################
65#ifndef CLP_MULTIPLE_FACTORIZATIONS   
66
67//-------------------------------------------------------------------
68// Default Constructor
69//-------------------------------------------------------------------
70ClpFactorization::ClpFactorization () :
71   CoinFactorization() 
72{
73#ifndef SLIM_CLP
74  networkBasis_ = NULL;
75#endif
76}
77
78//-------------------------------------------------------------------
79// Copy constructor
80//-------------------------------------------------------------------
81ClpFactorization::ClpFactorization (const ClpFactorization & rhs,
82                                    int dummyDenseIfSmaller) :
83   CoinFactorization(rhs) 
84{
85#ifndef SLIM_CLP
86  if (rhs.networkBasis_)
87    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
88  else
89    networkBasis_=NULL;
90#endif
91}
92
93ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
94   CoinFactorization(rhs) 
95{
96#ifndef SLIM_CLP
97  networkBasis_=NULL;
98#endif
99}
100
101//-------------------------------------------------------------------
102// Destructor
103//-------------------------------------------------------------------
104ClpFactorization::~ClpFactorization () 
105{
106#ifndef SLIM_CLP
107  delete networkBasis_;
108#endif
109}
110
111//----------------------------------------------------------------
112// Assignment operator
113//-------------------------------------------------------------------
114ClpFactorization &
115ClpFactorization::operator=(const ClpFactorization& rhs)
116{
117#ifdef CLP_FACTORIZATION_INSTRUMENT
118  factorization_instrument(-1);
119#endif
120  if (this != &rhs) {
121    CoinFactorization::operator=(rhs);
122#ifndef SLIM_CLP
123    delete networkBasis_;
124    if (rhs.networkBasis_)
125      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
126    else
127      networkBasis_=NULL;
128#endif
129  }
130#ifdef CLP_FACTORIZATION_INSTRUMENT
131  factorization_instrument(1);
132#endif
133  return *this;
134}
135#if 0
136static unsigned int saveList[10000];
137int numberSave=-1;
138inline bool isDense(int i) {
139  return ((saveList[i>>5]>>(i&31))&1)!=0;
140}
141inline void setDense(int i) {
142  unsigned int & value = saveList[i>>5];
143  int bit = i&31;
144  value |= (1<<bit);
145}
146#endif
147int 
148ClpFactorization::factorize ( ClpSimplex * model,
149                              int solveType, bool valuesPass)
150{
151#ifdef CLP_FACTORIZATION_INSTRUMENT
152  factorization_instrument(-1);
153#endif
154  ClpMatrixBase * matrix = model->clpMatrix(); 
155  int numberRows = model->numberRows();
156  int numberColumns = model->numberColumns();
157  if (!numberRows)
158    return 0;
159  // If too many compressions increase area
160  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
161    areaFactor_ *= 1.1;
162  }
163  //int numberPivots=numberPivots_;
164#if 0
165  if (model->algorithm()>0)
166    numberSave=-1;
167#endif
168#ifndef SLIM_CLP
169  if (!networkBasis_||doCheck) {
170#endif
171    status_=-99;
172    int * pivotVariable = model->pivotVariable();
173    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
174    while (status_<-98) {
175     
176      int i;
177      int numberBasic=0;
178      int numberRowBasic;
179      // Move pivot variables across if they look good
180      int * pivotTemp = model->rowArray(0)->getIndices();
181      assert (!model->rowArray(0)->getNumElements());
182      if (!matrix->rhsOffset(model)) {
183#if 0
184        if (numberSave>0) {
185          int nStill=0;
186          int nAtBound=0;
187          int nZeroDual=0;
188          CoinIndexedVector * array = model->rowArray(3);
189          CoinIndexedVector * objArray = model->columnArray(1);
190          array->clear();
191          objArray->clear();
192          double * cost = model->costRegion();
193          double tolerance = model->primalTolerance();
194          double offset=0.0;
195          for (i=0;i<numberRows;i++) {
196            int iPivot = pivotVariable[i];
197            if (iPivot<numberColumns&&isDense(iPivot)) {
198              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
199                nStill++;
200                double value=model->solutionRegion()[iPivot];
201                double dual = model->dualRowSolution()[i];
202                double lower=model->lowerRegion()[iPivot];
203                double upper=model->upperRegion()[iPivot];
204                ClpSimplex::Status status;
205                if (fabs(value-lower)<tolerance) {
206                  status=ClpSimplex::atLowerBound;
207                  nAtBound++;
208                } else if (fabs(value-upper)<tolerance) {
209                  nAtBound++;
210                  status=ClpSimplex::atUpperBound;
211                } else if (value>lower&&value<upper) {
212                  status=ClpSimplex::superBasic;
213                } else {
214                  status=ClpSimplex::basic;
215                }
216                if (status!=ClpSimplex::basic) {
217                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
218                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
219                    model->setRowStatus(i,ClpSimplex::basic);
220                    pivotVariable[i]=i+numberColumns;
221                    model->dualRowSolution()[i]=0.0;
222                    model->djRegion(0)[i]=0.0;
223                    array->add(i,dual);
224                    offset += dual*model->solutionRegion(0)[i];
225                  }
226                }
227                if (fabs(dual)<1.0e-5)
228                  nZeroDual++;
229              }
230            }
231          }
232          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
233                 numberSave,nStill,nAtBound,nZeroDual,offset);
234          if (array->getNumElements()) {
235            // modify costs
236            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
237                                               objArray);
238            array->clear();
239            int n=objArray->getNumElements();
240            int * indices = objArray->getIndices();
241            double * elements = objArray->denseVector();
242            for (i=0;i<n;i++) {
243              int iColumn = indices[i];
244              cost[iColumn] -= elements[iColumn];
245              elements[iColumn]=0.0;
246            }
247            objArray->setNumElements(0);
248          }
249        }
250#endif
251        // Seems to prefer things in order so quickest
252        // way is to go though like this
253        for (i=0;i<numberRows;i++) {
254          if (model->getRowStatus(i) == ClpSimplex::basic) 
255            pivotTemp[numberBasic++]=i;
256        }
257        numberRowBasic=numberBasic;
258        /* Put column basic variables into pivotVariable
259           This is done by ClpMatrixBase to allow override for gub
260        */
261        matrix->generalExpanded(model,0,numberBasic);
262      } else {
263        // Long matrix - do a different way
264        bool fullSearch=false;
265        for (i=0;i<numberRows;i++) {
266          int iPivot = pivotVariable[i];
267          if (iPivot>=numberColumns) {
268            pivotTemp[numberBasic++]=iPivot-numberColumns;
269          }
270        }
271        numberRowBasic=numberBasic;
272        for (i=0;i<numberRows;i++) {
273          int iPivot = pivotVariable[i];
274          if (iPivot<numberColumns) {
275            if (iPivot>=0) {
276              pivotTemp[numberBasic++]=iPivot;
277            } else {
278              // not full basis
279              fullSearch=true;
280              break;
281            }
282          }
283        }
284        if (fullSearch) {
285          // do slow way
286          numberBasic=0;
287          for (i=0;i<numberRows;i++) {
288            if (model->getRowStatus(i) == ClpSimplex::basic) 
289              pivotTemp[numberBasic++]=i;
290          }
291          numberRowBasic=numberBasic;
292          /* Put column basic variables into pivotVariable
293             This is done by ClpMatrixBase to allow override for gub
294          */
295          matrix->generalExpanded(model,0,numberBasic);
296        }
297      }
298      if (numberBasic>model->maximumBasic()) {
299#if 0 // ndef NDEBUG
300        printf("%d basic - should only be %d\n",
301               numberBasic,numberRows);
302#endif
303        // Take out some
304        numberBasic=numberRowBasic;
305        for (int i=0;i<numberColumns;i++) {
306          if (model->getColumnStatus(i) == ClpSimplex::basic) {
307            if (numberBasic<numberRows)
308              numberBasic++;
309            else
310              model->setColumnStatus(i,ClpSimplex::superBasic);
311          }
312        }
313        numberBasic=numberRowBasic;
314        matrix->generalExpanded(model,0,numberBasic);
315      }
316#ifndef SLIM_CLP
317      // see if matrix a network
318#ifndef NO_RTTI
319      ClpNetworkMatrix* networkMatrix =
320        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
321#else
322      ClpNetworkMatrix* networkMatrix = NULL;
323      if (model->clpMatrix()->type()==11)
324        networkMatrix = 
325        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
326#endif
327      // If network - still allow ordinary factorization first time for laziness
328      if (networkMatrix)
329        biasLU_=0; // All to U if network
330      //int saveMaximumPivots = maximumPivots();
331      delete networkBasis_;
332      networkBasis_ = NULL;
333      if (networkMatrix&&!doCheck)
334        maximumPivots(1);
335#endif
336      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
337      while (status_==-99) {
338        // maybe for speed will be better to leave as many regions as possible
339        gutsOfDestructor();
340        gutsOfInitialize(2);
341        CoinBigIndex numberElements=numberRowBasic;
342
343        // compute how much in basis
344
345        int i;
346        // can change for gub
347        int numberColumnBasic = numberBasic-numberRowBasic;
348
349        numberElements +=matrix->countBasis(model,
350                                           pivotTemp+numberRowBasic, 
351                                           numberRowBasic,
352                                            numberColumnBasic);
353        // and recompute as network side say different
354        if (model->numberIterations())
355          numberRowBasic = numberBasic - numberColumnBasic;
356        numberElements = 3 * numberBasic + 3 * numberElements + 20000;
357#if 0
358        // If iteration not zero then may be compressed
359        getAreas ( !model->numberIterations() ? numberRows : numberBasic,
360                   numberRowBasic+numberColumnBasic, numberElements,
361                   2 * numberElements );
362#else
363        getAreas ( numberRows,
364                   numberRowBasic+numberColumnBasic, numberElements,
365                   2 * numberElements );
366#endif
367        //fill
368        // Fill in counts so we can skip part of preProcess
369        int * numberInRow = numberInRow_.array();
370        int * numberInColumn = numberInColumn_.array();
371        CoinZeroN ( numberInRow, numberRows_ + 1 );
372        CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 );
373        double * elementU = elementU_.array();
374        int * indexRowU = indexRowU_.array();
375        CoinBigIndex * startColumnU = startColumnU_.array();
376        for (i=0;i<numberRowBasic;i++) {
377          int iRow = pivotTemp[i];
378          indexRowU[i]=iRow;
379          startColumnU[i]=i;
380          elementU[i]=slackValue_;
381          numberInRow[iRow]=1;
382          numberInColumn[i]=1;
383        }
384        startColumnU[numberRowBasic]=numberRowBasic;
385        // can change for gub so redo
386        numberColumnBasic = numberBasic-numberRowBasic;
387        matrix->fillBasis(model, 
388                          pivotTemp+numberRowBasic, 
389                          numberColumnBasic,
390                          indexRowU, 
391                          startColumnU+numberRowBasic,
392                          numberInRow,
393                          numberInColumn+numberRowBasic,
394                          elementU);
395#if 0
396        {
397          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
398          for (int i=0;i<numberElements;i++)
399            printf("row %d col %d value %g\n",indexRowU_.array()[i],indexColumnU_[i],
400                   elementU_.array()[i]);
401        }
402#endif
403        // recompute number basic
404        numberBasic = numberRowBasic+numberColumnBasic;
405        if (numberBasic) 
406          numberElements = startColumnU[numberBasic-1]
407            +numberInColumn[numberBasic-1];
408        else
409          numberElements=0;
410        lengthU_ = numberElements;
411        //saveFactorization("dump.d");
412        if (biasLU_>=3||numberRows_!=numberColumns_)
413          preProcess ( 2 );
414        else
415          preProcess ( 3 ); // no row copy
416        factor (  );
417        if (status_==-99) {
418          // get more memory
419          areaFactor(2.0*areaFactor());
420        } else if (status_==-1&&model->numberIterations()==0&&
421                   denseThreshold_) {
422          // Round again without dense
423          denseThreshold_=0;
424          status_=-99;
425        }
426      }
427      // If we get here status is 0 or -1
428      if (status_ == 0) {
429        // We may need to tamper with order and redo - e.g. network with side
430        int useNumberRows = numberRows;
431        // **** we will also need to add test in dual steepest to do
432        // as we do for network
433        matrix->generalExpanded(model,12,useNumberRows);
434        const int * permuteBack = permuteBack_.array();
435        const int * back = pivotColumnBack_.array();
436        //int * pivotTemp = pivotColumn_.array();
437        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
438        // Redo pivot order
439        for (i=0;i<numberRowBasic;i++) {
440          int k = pivotTemp[i];
441          // so rowIsBasic[k] would be permuteBack[back[i]]
442          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
443        }
444        for (;i<useNumberRows;i++) {
445          int k = pivotTemp[i];
446          // so rowIsBasic[k] would be permuteBack[back[i]]
447          pivotVariable[permuteBack[back[i]]]=k;
448        }
449#if 0
450        if (numberSave>=0) {
451          numberSave=numberDense_;
452          memset(saveList,0,((numberRows_+31)>>5)*sizeof(int));
453          for (i=numberRows_-numberSave;i<numberRows_;i++) {
454            int k=pivotTemp[pivotColumn_.array()[i]];
455            setDense(k);
456          }
457        }
458#endif
459        // Set up permutation vector
460        // these arrays start off as copies of permute
461        // (and we could use permute_ instead of pivotColumn (not back though))
462        ClpDisjointCopyN ( permute_.array(), useNumberRows , pivotColumn_.array()  );
463        ClpDisjointCopyN ( permuteBack_.array(), useNumberRows , pivotColumnBack_.array()  );
464#ifndef SLIM_CLP
465        if (networkMatrix) {
466          maximumPivots(CoinMax(2000,maximumPivots()));
467          // redo arrays
468          for (int iRow=0;iRow<4;iRow++) {
469            int length =model->numberRows()+maximumPivots();
470            if (iRow==3||model->objectiveAsObject()->type()>1)
471              length += model->numberColumns();
472            model->rowArray(iRow)->reserve(length);
473          }
474          // create network factorization
475          if (doCheck)
476            delete networkBasis_; // temp
477          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
478                                              pivotRegion_.array(),
479                                              permuteBack_.array(),
480                                              startColumnU_.array(),
481                                              numberInColumn_.array(),
482                                              indexRowU_.array(),
483                                              elementU_.array());
484          // kill off arrays in ordinary factorization
485          if (!doCheck) {
486            gutsOfDestructor();
487            // but make sure numberRows_ set
488            numberRows_ = model->numberRows();
489            status_=0;
490#if 0
491            // but put back permute arrays so odd things will work
492            int numberRows = model->numberRows();
493            pivotColumnBack_ = new int [numberRows];
494            permute_ = new int [numberRows];
495            int i;
496            for (i=0;i<numberRows;i++) {
497              pivotColumnBack_[i]=i;
498              permute_[i]=i;
499            }
500#endif
501          }
502        } else {
503#endif
504          // See if worth going sparse and when
505          if (numberFtranCounts_>100) {
506            ftranCountInput_= CoinMax(ftranCountInput_,1.0);
507            ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,1.0);
508            ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,1.0);
509            ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,1.0);
510            if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
511              btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,1.0);
512              btranAverageAfterR_ = CoinMax(btranCountAfterR_/btranCountAfterU_,1.0);
513              btranAverageAfterL_ = CoinMax(btranCountAfterL_/btranCountAfterR_,1.0);
514            } else {
515              // we have not done any useful btrans (values pass?)
516              btranAverageAfterU_ = 1.0;
517              btranAverageAfterR_ = 1.0;
518              btranAverageAfterL_ = 1.0;
519            }
520          }
521          // scale back
522         
523          ftranCountInput_ *= 0.8;
524          ftranCountAfterL_ *= 0.8;
525          ftranCountAfterR_ *= 0.8;
526          ftranCountAfterU_ *= 0.8;
527          btranCountInput_ *= 0.8;
528          btranCountAfterU_ *= 0.8;
529          btranCountAfterR_ *= 0.8;
530          btranCountAfterL_ *= 0.8;
531#ifndef SLIM_CLP
532        }
533#endif
534      } else if (status_ == -1&&(solveType==0||solveType==2)) {
535        // This needs redoing as it was merged coding - does not need array
536        int numberTotal = numberRows+numberColumns;
537        int * isBasic = new int [numberTotal];
538        int * rowIsBasic = isBasic+numberColumns;
539        int * columnIsBasic = isBasic;
540        for (i=0;i<numberTotal;i++) 
541          isBasic[i]=-1;
542        for (i=0;i<numberRowBasic;i++) {
543          int iRow = pivotTemp[i];
544          rowIsBasic[iRow]=1;
545        }
546        for (;i<numberBasic;i++) {
547          int iColumn = pivotTemp[i];
548          columnIsBasic[iColumn]=1;
549        }
550        numberBasic=0;
551        for (i=0;i<numberRows;i++) 
552          pivotVariable[i]=-1;
553        // mark as basic or non basic
554        const int * pivotColumn = pivotColumn_.array();
555        for (i=0;i<numberRows;i++) {
556          if (rowIsBasic[i]>=0) {
557            if (pivotColumn[numberBasic]>=0) {
558              rowIsBasic[i]=pivotColumn[numberBasic];
559            } else {
560              rowIsBasic[i]=-1;
561              model->setRowStatus(i,ClpSimplex::superBasic);
562            }
563            numberBasic++;
564          }
565        }
566        for (i=0;i<numberColumns;i++) {
567          if (columnIsBasic[i]>=0) {
568            if (pivotColumn[numberBasic]>=0) 
569              columnIsBasic[i]=pivotColumn[numberBasic];
570            else
571              columnIsBasic[i]=-1;
572            numberBasic++;
573          }
574        }
575        // leave pivotVariable in useful form for cleaning basis
576        int * pivotVariable = model->pivotVariable();
577        for (i=0;i<numberRows;i++) {
578          pivotVariable[i]=-1;
579        }
580
581        for (i=0;i<numberRows;i++) {
582          if (model->getRowStatus(i) == ClpSimplex::basic) {
583            int iPivot = rowIsBasic[i];
584            if (iPivot>=0) 
585              pivotVariable[iPivot]=i+numberColumns;
586          }
587        }
588        for (i=0;i<numberColumns;i++) {
589          if (model->getColumnStatus(i) == ClpSimplex::basic) {
590            int iPivot = columnIsBasic[i];
591            if (iPivot>=0) 
592              pivotVariable[iPivot]=i;
593          }
594        }
595        delete [] isBasic;
596        double * columnLower = model->lowerRegion();
597        double * columnUpper = model->upperRegion();
598        double * columnActivity = model->solutionRegion();
599        double * rowLower = model->lowerRegion(0);
600        double * rowUpper = model->upperRegion(0);
601        double * rowActivity = model->solutionRegion(0);
602        //redo basis - first take ALL columns out
603        int iColumn;
604        double largeValue = model->largeValue();
605        for (iColumn=0;iColumn<numberColumns;iColumn++) {
606          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
607            // take out
608            if (!valuesPass) {
609              double lower=columnLower[iColumn];
610              double upper=columnUpper[iColumn];
611              double value=columnActivity[iColumn];
612              if (lower>-largeValue||upper<largeValue) {
613                if (fabs(value-lower)<fabs(value-upper)) {
614                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
615                  columnActivity[iColumn]=lower;
616                } else {
617                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
618                  columnActivity[iColumn]=upper;
619                }
620              } else {
621                model->setColumnStatus(iColumn,ClpSimplex::isFree);
622              }
623            } else {
624              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
625            }
626          }
627        }
628        int iRow;
629        for (iRow=0;iRow<numberRows;iRow++) {
630          int iSequence=pivotVariable[iRow];
631          if (iSequence>=0) {
632            // basic
633            if (iSequence>=numberColumns) {
634              // slack in - leave
635              //assert (iSequence-numberColumns==iRow);
636            } else {
637              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
638              // put back structural
639              model->setColumnStatus(iSequence,ClpSimplex::basic);
640            }
641          } else {
642            // put in slack
643            model->setRowStatus(iRow,ClpSimplex::basic);
644          }
645        }
646        // Put back any key variables for gub (status_ not touched)
647        matrix->generalExpanded(model,1,status_);
648        // signal repeat
649        status_=-99;
650        // set fixed if they are
651        for (iRow=0;iRow<numberRows;iRow++) {
652          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
653            if (rowLower[iRow]==rowUpper[iRow]) {
654              rowActivity[iRow]=rowLower[iRow];
655              model->setRowStatus(iRow,ClpSimplex::isFixed);
656            }
657          }
658        }
659        for (iColumn=0;iColumn<numberColumns;iColumn++) {
660          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
661            if (columnLower[iColumn]==columnUpper[iColumn]) {
662              columnActivity[iColumn]=columnLower[iColumn];
663              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
664            }
665          }
666        }
667      } 
668    }
669#ifndef SLIM_CLP
670  } else {
671    // network - fake factorization - do nothing
672    status_=0;
673    numberPivots_ = 0;
674  }
675#endif
676#ifndef SLIM_CLP
677  if (!status_) {
678    // take out part if quadratic
679    if (model->algorithm()==2) {
680      ClpObjective * obj = model->objectiveAsObject();
681#ifndef NDEBUG
682      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
683      assert (quadraticObj);
684#else
685      ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
686#endif
687      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
688      int numberXColumns = quadratic->getNumCols();
689      assert (numberXColumns<numberColumns);
690      int base = numberColumns-numberXColumns;
691      int * which = new int [numberXColumns];
692      int * pivotVariable = model->pivotVariable();
693      int * permute = pivotColumn();
694      int i;
695      int n=0;
696      for (i=0;i<numberRows;i++) {
697        int iSj = pivotVariable[i]-base;
698        if (iSj>=0&&iSj<numberXColumns) 
699          which[n++]=permute[i];
700      }
701      if (n)
702        emptyRows(n,which);
703      delete [] which;
704    }
705  }
706#endif
707#ifdef CLP_FACTORIZATION_INSTRUMENT
708  factorization_instrument(2);
709#endif
710  return status_;
711}
712/* Replaces one Column to basis,
713   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
714   If checkBeforeModifying is true will do all accuracy checks
715   before modifying factorization.  Whether to set this depends on
716   speed considerations.  You could just do this on first iteration
717   after factorization and thereafter re-factorize
718   partial update already in U */
719int 
720ClpFactorization::replaceColumn ( const ClpSimplex * model, 
721                                  CoinIndexedVector * regionSparse,
722                                  CoinIndexedVector * tableauColumn,
723                                  int pivotRow,
724                                  double pivotCheck ,
725                                  bool checkBeforeModifying)
726{
727  int returnCode;
728#ifndef SLIM_CLP
729  if (!networkBasis_) {
730#endif
731#ifdef CLP_FACTORIZATION_INSTRUMENT
732    factorization_instrument(-1);
733#endif
734    // see if FT
735    if (doForrestTomlin_) {
736      returnCode= CoinFactorization::replaceColumn(regionSparse,
737                                              pivotRow,
738                                              pivotCheck,
739                                              checkBeforeModifying);
740    } else {
741      returnCode= CoinFactorization::replaceColumnPFI(tableauColumn,
742                                              pivotRow,pivotCheck); // Note array
743    }
744#ifdef CLP_FACTORIZATION_INSTRUMENT
745    factorization_instrument(3);
746#endif
747
748#ifndef SLIM_CLP
749  } else {
750    if (doCheck) {
751      returnCode = CoinFactorization::replaceColumn(regionSparse,
752                                                        pivotRow,
753                                                        pivotCheck,
754                                                        checkBeforeModifying);
755      networkBasis_->replaceColumn(regionSparse,
756                                   pivotRow);
757    } else {
758      // increase number of pivots
759      numberPivots_++;
760      returnCode = networkBasis_->replaceColumn(regionSparse,
761                                   pivotRow);
762    }
763  }
764#endif
765  return returnCode;
766}
767
768/* Updates one column (FTRAN) from region2
769   number returned is negative if no room
770   region1 starts as zero and is zero at end */
771int 
772ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
773                                   CoinIndexedVector * regionSparse2)
774{
775#ifdef CLP_DEBUG
776  regionSparse->checkClear();
777#endif
778  if (!numberRows_)
779    return 0;
780#ifndef SLIM_CLP
781  if (!networkBasis_) {
782#endif
783#ifdef CLP_FACTORIZATION_INSTRUMENT
784    factorization_instrument(-1);
785#endif
786    collectStatistics_ = true;
787    int returnCode= CoinFactorization::updateColumnFT(regionSparse,
788                                                       regionSparse2);
789    collectStatistics_ = false;
790#ifdef CLP_FACTORIZATION_INSTRUMENT
791    factorization_instrument(4);
792#endif
793    return returnCode;
794#ifndef SLIM_CLP
795  } else {
796#ifdef CHECK_NETWORK
797    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
798    double * check = new double[numberRows_];
799    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
800                                                       regionSparse2);
801    networkBasis_->updateColumn(regionSparse,save,-1);
802    int i;
803    double * array = regionSparse2->denseVector();
804    int * indices = regionSparse2->getIndices();
805    int n=regionSparse2->getNumElements();
806    memset(check,0,numberRows_*sizeof(double));
807    double * array2 = save->denseVector();
808    int * indices2 = save->getIndices();
809    int n2=save->getNumElements();
810    assert (n==n2);
811    if (save->packedMode()) {
812      for (i=0;i<n;i++) {
813        check[indices[i]]=array[i];
814      }
815      for (i=0;i<n;i++) {
816        double value2 = array2[i];
817        assert (check[indices2[i]]==value2);
818      }
819    } else {
820      for (i=0;i<numberRows_;i++) {
821        double value1 = array[i];
822        double value2 = array2[i];
823        assert (value1==value2);
824      }
825    }
826    delete save;
827    delete [] check;
828    return returnCode;
829#else
830    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
831    return 1;
832#endif
833  }
834#endif
835}
836/* Updates one column (FTRAN) from region2
837   number returned is negative if no room
838   region1 starts as zero and is zero at end */
839int 
840ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
841                                 CoinIndexedVector * regionSparse2,
842                                 bool noPermute) const
843{
844#ifdef CLP_DEBUG
845  if (!noPermute)
846    regionSparse->checkClear();
847#endif
848  if (!numberRows_)
849    return 0;
850#ifndef SLIM_CLP
851  if (!networkBasis_) {
852#endif
853#ifdef CLP_FACTORIZATION_INSTRUMENT
854    factorization_instrument(-1);
855#endif
856    collectStatistics_ = true;
857    int returnCode= CoinFactorization::updateColumn(regionSparse,
858                                                     regionSparse2,
859                                                     noPermute);
860    collectStatistics_ = false;
861#ifdef CLP_FACTORIZATION_INSTRUMENT
862    factorization_instrument(5);
863#endif
864    return returnCode;
865#ifndef SLIM_CLP
866  } else {
867#ifdef CHECK_NETWORK
868    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
869    double * check = new double[numberRows_];
870    int returnCode = CoinFactorization::updateColumn(regionSparse,
871                                                     regionSparse2,
872                                                     noPermute);
873    networkBasis_->updateColumn(regionSparse,save,-1);
874    int i;
875    double * array = regionSparse2->denseVector();
876    int * indices = regionSparse2->getIndices();
877    int n=regionSparse2->getNumElements();
878    memset(check,0,numberRows_*sizeof(double));
879    double * array2 = save->denseVector();
880    int * indices2 = save->getIndices();
881    int n2=save->getNumElements();
882    assert (n==n2);
883    if (save->packedMode()) {
884      for (i=0;i<n;i++) {
885        check[indices[i]]=array[i];
886      }
887      for (i=0;i<n;i++) {
888        double value2 = array2[i];
889        assert (check[indices2[i]]==value2);
890      }
891    } else {
892      for (i=0;i<numberRows_;i++) {
893        double value1 = array[i];
894        double value2 = array2[i];
895        assert (value1==value2);
896      }
897    }
898    delete save;
899    delete [] check;
900    return returnCode;
901#else
902    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
903    return 1;
904#endif
905  }
906#endif
907}
908/* Updates one column (FTRAN) from region2
909   Tries to do FT update
910   number returned is negative if no room.
911   Also updates region3
912   region1 starts as zero and is zero at end */
913int 
914ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
915                                       CoinIndexedVector * regionSparse2,
916                                       CoinIndexedVector * regionSparse3,
917                                       bool noPermuteRegion3) 
918{
919  int returnCode=updateColumnFT(regionSparse1,regionSparse2);
920  updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
921  return returnCode;
922}
923/* Updates one column (FTRAN) from region2
924   number returned is negative if no room
925   region1 starts as zero and is zero at end */
926int 
927ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
928                                 CoinIndexedVector * regionSparse2,
929                                 bool noPermute) const
930{
931  if (!noPermute)
932    regionSparse->checkClear();
933  if (!numberRows_)
934    return 0;
935  collectStatistics_ = false;
936  int returnCode= CoinFactorization::updateColumn(regionSparse,
937                                                   regionSparse2,
938                                                   noPermute);
939  return returnCode;
940}
941/* Updates one column (BTRAN) from region2
942   region1 starts as zero and is zero at end */
943int 
944ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
945                          CoinIndexedVector * regionSparse2) const
946{
947  if (!numberRows_)
948    return 0;
949#ifndef SLIM_CLP
950  if (!networkBasis_) {
951#endif
952#ifdef CLP_FACTORIZATION_INSTRUMENT
953    factorization_instrument(-1);
954#endif
955    collectStatistics_ = true;
956    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
957                                                    regionSparse2);
958    collectStatistics_ = false;
959#ifdef CLP_FACTORIZATION_INSTRUMENT
960    factorization_instrument(6);
961#endif
962    return returnCode;
963#ifndef SLIM_CLP
964  } else {
965#ifdef CHECK_NETWORK
966    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
967    double * check = new double[numberRows_];
968    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
969                                                              regionSparse2);
970    networkBasis_->updateColumnTranspose(regionSparse,save);
971    int i;
972    double * array = regionSparse2->denseVector();
973    int * indices = regionSparse2->getIndices();
974    int n=regionSparse2->getNumElements();
975    memset(check,0,numberRows_*sizeof(double));
976    double * array2 = save->denseVector();
977    int * indices2 = save->getIndices();
978    int n2=save->getNumElements();
979    assert (n==n2);
980    if (save->packedMode()) {
981      for (i=0;i<n;i++) {
982        check[indices[i]]=array[i];
983      }
984      for (i=0;i<n;i++) {
985        double value2 = array2[i];
986        assert (check[indices2[i]]==value2);
987      }
988    } else {
989      for (i=0;i<numberRows_;i++) {
990        double value1 = array[i];
991        double value2 = array2[i];
992        assert (value1==value2);
993      }
994    }
995    delete save;
996    delete [] check;
997    return returnCode;
998#else
999    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
1000#endif
1001  }
1002#endif
1003}
1004/* makes a row copy of L for speed and to allow very sparse problems */
1005void 
1006ClpFactorization::goSparse()
1007{
1008#ifdef CLP_FACTORIZATION_INSTRUMENT
1009  factorization_instrument(-1);
1010#endif
1011#ifndef SLIM_CLP
1012  if (!networkBasis_) 
1013#endif
1014    CoinFactorization::goSparse();
1015#ifdef CLP_FACTORIZATION_INSTRUMENT
1016  factorization_instrument(7);
1017#endif
1018}
1019// Cleans up i.e. gets rid of network basis
1020void 
1021ClpFactorization::cleanUp()
1022{
1023#ifndef SLIM_CLP
1024  delete networkBasis_;
1025  networkBasis_=NULL;
1026#endif
1027  resetStatistics();
1028}
1029/// Says whether to redo pivot order
1030bool 
1031ClpFactorization::needToReorder() const
1032{
1033#ifdef CHECK_NETWORK
1034  return true;
1035#endif
1036#ifndef SLIM_CLP
1037  if (!networkBasis_)
1038#endif
1039    return true;
1040#ifndef SLIM_CLP
1041  else
1042    return false;
1043#endif
1044}
1045// Get weighted row list
1046void
1047ClpFactorization::getWeights(int * weights) const
1048{
1049#ifdef CLP_FACTORIZATION_INSTRUMENT
1050  factorization_instrument(-1);
1051#endif
1052#ifndef SLIM_CLP
1053  if (networkBasis_) {
1054    // Network - just unit
1055    for (int i=0;i<numberRows_;i++) 
1056      weights[i]=1;
1057    return;
1058  }
1059#endif
1060  int * numberInRow = numberInRow_.array();
1061  int * numberInColumn = numberInColumn_.array();
1062  int * permuteBack = pivotColumnBack_.array();
1063  int * indexRowU = indexRowU_.array();
1064  const CoinBigIndex * startColumnU = startColumnU_.array();
1065  const CoinBigIndex * startRowL = startRowL_.array();
1066  if (!startRowL||!numberInRow_.array()) {
1067    int * temp = new int[numberRows_];
1068    memset(temp,0,numberRows_*sizeof(int));
1069    int i;
1070    for (i=0;i<numberRows_;i++) {
1071      // one for pivot
1072      temp[i]++;
1073      CoinBigIndex j;
1074      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
1075        int iRow=indexRowU[j];
1076        temp[iRow]++;
1077      }
1078    }
1079    CoinBigIndex * startColumnL = startColumnL_.array();
1080    int * indexRowL = indexRowL_.array();
1081    for (i=baseL_;i<baseL_+numberL_;i++) {
1082      CoinBigIndex j;
1083      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
1084        int iRow = indexRowL[j];
1085        temp[iRow]++;
1086      }
1087    }
1088    for (i=0;i<numberRows_;i++) {
1089      int number = temp[i];
1090      int iPermute = permuteBack[i];
1091      weights[iPermute]=number;
1092    }
1093    delete [] temp;
1094  } else {
1095    int i;
1096    for (i=0;i<numberRows_;i++) {
1097      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
1098      //number = startRowL[i+1]-startRowL[i]+1;
1099      //number = numberInRow[i]+1;
1100      int iPermute = permuteBack[i];
1101      weights[iPermute]=number;
1102    }
1103  }
1104#ifdef CLP_FACTORIZATION_INSTRUMENT
1105  factorization_instrument(8);
1106#endif
1107}
1108#else
1109// This one allows multiple factorizations
1110
1111//-------------------------------------------------------------------
1112// Default Constructor
1113//-------------------------------------------------------------------
1114ClpFactorization::ClpFactorization () 
1115{
1116#ifndef SLIM_CLP
1117  networkBasis_ = NULL;
1118#endif
1119  //coinFactorizationA_ = NULL;
1120  coinFactorizationA_ = new CoinFactorization() ;
1121  coinFactorizationB_ = NULL;
1122  //coinFactorizationB_ = new CoinSmallFactorization();
1123  forceB_=0;
1124  goOslThreshold_ = -1;
1125  goDenseThreshold_ = -1;
1126  goSmallThreshold_ = -1;
1127}
1128
1129//-------------------------------------------------------------------
1130// Copy constructor
1131//-------------------------------------------------------------------
1132ClpFactorization::ClpFactorization (const ClpFactorization & rhs,
1133                                    int denseIfSmaller) 
1134{
1135#ifdef CLP_FACTORIZATION_INSTRUMENT
1136  factorization_instrument(-1);
1137#endif
1138#ifndef SLIM_CLP
1139  if (rhs.networkBasis_)
1140    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1141  else
1142    networkBasis_=NULL;
1143#endif
1144  forceB_ = rhs.forceB_;
1145  goOslThreshold_ = rhs.goOslThreshold_;
1146  goDenseThreshold_ = rhs.goDenseThreshold_;
1147  goSmallThreshold_ = rhs.goSmallThreshold_;
1148  int goDense = 0;
1149  if (denseIfSmaller>0&&!rhs.coinFactorizationB_) {
1150    if (denseIfSmaller<=goDenseThreshold_) 
1151      goDense=1;
1152    else if (denseIfSmaller<=goSmallThreshold_) 
1153      goDense=2;
1154    else if (denseIfSmaller<=goOslThreshold_) 
1155      goDense=3;
1156  } else if (denseIfSmaller<0) {
1157    if (-denseIfSmaller<=goDenseThreshold_) 
1158      goDense=1;
1159    else if (-denseIfSmaller<=goSmallThreshold_) 
1160      goDense=2;
1161    else if (-denseIfSmaller<=goOslThreshold_) 
1162      goDense=3;
1163  }
1164  if (rhs.coinFactorizationA_&&!goDense)
1165    coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
1166  else
1167    coinFactorizationA_=NULL;
1168  if (rhs.coinFactorizationB_&&(denseIfSmaller>=0||!goDense))
1169    coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1170  else
1171    coinFactorizationB_=NULL;
1172  if (goDense) {
1173    delete coinFactorizationB_;
1174    if (goDense==1)
1175      coinFactorizationB_ = new CoinDenseFactorization();
1176    else if (goDense==2)
1177      coinFactorizationB_ = new CoinSimpFactorization();
1178    else
1179      coinFactorizationB_ = new CoinOslFactorization();
1180    if (rhs.coinFactorizationA_) {
1181      coinFactorizationB_->maximumPivots(rhs.coinFactorizationA_->maximumPivots());
1182      coinFactorizationB_->pivotTolerance(rhs.coinFactorizationA_->pivotTolerance());
1183      coinFactorizationB_->zeroTolerance(rhs.coinFactorizationA_->zeroTolerance());
1184    } else {
1185      assert (coinFactorizationB_);
1186      coinFactorizationB_->maximumPivots(rhs.coinFactorizationB_->maximumPivots());
1187      coinFactorizationB_->pivotTolerance(rhs.coinFactorizationB_->pivotTolerance());
1188      coinFactorizationB_->zeroTolerance(rhs.coinFactorizationB_->zeroTolerance());
1189    }
1190  }
1191  assert (!coinFactorizationA_||!coinFactorizationB_);
1192#ifdef CLP_FACTORIZATION_INSTRUMENT
1193  factorization_instrument(1);
1194#endif
1195}
1196
1197ClpFactorization::ClpFactorization (const CoinFactorization & rhs) 
1198{
1199#ifdef CLP_FACTORIZATION_INSTRUMENT
1200  factorization_instrument(-1);
1201#endif
1202#ifndef SLIM_CLP
1203  networkBasis_=NULL;
1204#endif
1205  coinFactorizationA_ = new CoinFactorization(rhs);
1206  coinFactorizationB_=NULL;
1207#ifdef CLP_FACTORIZATION_INSTRUMENT
1208  factorization_instrument(1);
1209#endif
1210  forceB_=0;
1211  goOslThreshold_ = -1;
1212  goDenseThreshold_ = -1;
1213  goSmallThreshold_ = -1;
1214  assert (!coinFactorizationA_||!coinFactorizationB_);
1215}
1216
1217ClpFactorization::ClpFactorization (const CoinSmallFactorization & rhs) 
1218{
1219#ifdef CLP_FACTORIZATION_INSTRUMENT
1220  factorization_instrument(-1);
1221#endif
1222#ifndef SLIM_CLP
1223  networkBasis_=NULL;
1224#endif
1225  coinFactorizationA_ = NULL;
1226  coinFactorizationB_ = rhs.clone();
1227  //coinFactorizationB_ = new CoinSmallFactorization(rhs);
1228  forceB_=0;
1229  goOslThreshold_ = -1;
1230  goDenseThreshold_ = -1;
1231  goSmallThreshold_ = -1;
1232#ifdef CLP_FACTORIZATION_INSTRUMENT
1233  factorization_instrument(1);
1234#endif
1235  assert (!coinFactorizationA_||!coinFactorizationB_);
1236}
1237
1238//-------------------------------------------------------------------
1239// Destructor
1240//-------------------------------------------------------------------
1241ClpFactorization::~ClpFactorization () 
1242{
1243#ifndef SLIM_CLP
1244  delete networkBasis_;
1245#endif
1246  delete coinFactorizationA_;
1247  delete coinFactorizationB_;
1248}
1249
1250//----------------------------------------------------------------
1251// Assignment operator
1252//-------------------------------------------------------------------
1253ClpFactorization &
1254ClpFactorization::operator=(const ClpFactorization& rhs)
1255{
1256#ifdef CLP_FACTORIZATION_INSTRUMENT
1257  factorization_instrument(-1);
1258#endif
1259  if (this != &rhs) {
1260#ifndef SLIM_CLP
1261    delete networkBasis_;
1262    if (rhs.networkBasis_)
1263      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1264    else
1265      networkBasis_=NULL;
1266#endif
1267    forceB_ = rhs.forceB_;
1268    goOslThreshold_ = rhs.goOslThreshold_;
1269    goDenseThreshold_ = rhs.goDenseThreshold_;
1270    goSmallThreshold_ = rhs.goSmallThreshold_;
1271    if (rhs.coinFactorizationA_) {
1272      if (coinFactorizationA_)
1273        *coinFactorizationA_ = *(rhs.coinFactorizationA_);
1274      else
1275        coinFactorizationA_ = new CoinFactorization(*rhs.coinFactorizationA_);
1276    } else {
1277      delete coinFactorizationA_;
1278      coinFactorizationA_=NULL;
1279    }
1280   
1281    if (rhs.coinFactorizationB_) {
1282      if (coinFactorizationB_) {
1283        CoinDenseFactorization * denseR = dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
1284        CoinDenseFactorization * dense = dynamic_cast<CoinDenseFactorization *>(coinFactorizationB_);
1285        CoinOslFactorization * oslR = dynamic_cast<CoinOslFactorization *>(rhs.coinFactorizationB_);
1286        CoinOslFactorization * osl = dynamic_cast<CoinOslFactorization *>(coinFactorizationB_);
1287        CoinSimpFactorization * simpR = dynamic_cast<CoinSimpFactorization *>(rhs.coinFactorizationB_);
1288        CoinSimpFactorization * simp = dynamic_cast<CoinSimpFactorization *>(coinFactorizationB_);
1289        if (dense&&denseR) {
1290          *dense=*denseR;
1291        } else if (osl&&oslR) {
1292          *osl=*oslR;
1293        } else if (simp&&simpR) {
1294          *simp=*simpR;
1295        } else {
1296          delete coinFactorizationB_;
1297          coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1298        }
1299      } else {
1300        coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1301      }
1302    } else {
1303      delete coinFactorizationB_;
1304      coinFactorizationB_=NULL;
1305    }
1306  }
1307#ifdef CLP_FACTORIZATION_INSTRUMENT
1308  factorization_instrument(1);
1309#endif
1310  assert (!coinFactorizationA_||!coinFactorizationB_);
1311  return *this;
1312}
1313// Go over to dense code
1314void 
1315ClpFactorization::goDenseOrSmall(int numberRows) 
1316{
1317  if (!forceB_) {
1318    if (numberRows<=goDenseThreshold_) {
1319      delete coinFactorizationA_;
1320      delete coinFactorizationB_;
1321      coinFactorizationA_ = NULL;
1322      coinFactorizationB_ = new CoinDenseFactorization();
1323      //printf("going dense\n");
1324    } else if (numberRows<=goSmallThreshold_) {
1325      delete coinFactorizationA_;
1326      delete coinFactorizationB_;
1327      coinFactorizationA_ = NULL;
1328      coinFactorizationB_ = new CoinSimpFactorization();
1329      //printf("going small\n");
1330    } else if (numberRows<=goOslThreshold_) {
1331      delete coinFactorizationA_;
1332      delete coinFactorizationB_;
1333      coinFactorizationA_ = NULL;
1334      coinFactorizationB_ = new CoinOslFactorization();
1335      //printf("going small\n");
1336    }
1337  }
1338  assert (!coinFactorizationA_||!coinFactorizationB_);
1339}
1340// If nonzero force use of 1,dense 2,small 3,osl
1341void 
1342ClpFactorization::forceOtherFactorization(int which)
1343{
1344  delete coinFactorizationB_;
1345  forceB_=0;
1346  coinFactorizationB_ = NULL;
1347  if (which>0&&which<4) {
1348    delete coinFactorizationA_;
1349    coinFactorizationA_ = NULL;
1350    forceB_=which;
1351    switch (which) {
1352    case 1:
1353      coinFactorizationB_ = new CoinDenseFactorization();
1354      goDenseThreshold_=COIN_INT_MAX;
1355    break;
1356    case 2:
1357      coinFactorizationB_ = new CoinSimpFactorization();
1358      goSmallThreshold_=COIN_INT_MAX;
1359    break;
1360    case 3:
1361      coinFactorizationB_ = new CoinOslFactorization();
1362      goOslThreshold_=COIN_INT_MAX;
1363    break;
1364    }
1365  } else if (!coinFactorizationA_) {
1366    coinFactorizationA_ = new CoinFactorization();
1367  }
1368}
1369int 
1370ClpFactorization::factorize ( ClpSimplex * model,
1371                              int solveType, bool valuesPass)
1372{
1373  //if ((model->specialOptions()&16384))
1374  //printf("factor at %d iterations\n",model->numberIterations());
1375  ClpMatrixBase * matrix = model->clpMatrix(); 
1376  int numberRows = model->numberRows();
1377  int numberColumns = model->numberColumns();
1378  if (!numberRows)
1379    return 0;
1380#ifdef CLP_FACTORIZATION_INSTRUMENT
1381  factorization_instrument(-1);
1382#endif
1383  bool anyChanged=false;
1384  if (coinFactorizationB_) {
1385    coinFactorizationB_->setStatus(-99);
1386    int * pivotVariable = model->pivotVariable();
1387    //returns 0 -okay, -1 singular, -2 too many in basis */
1388    // allow dense
1389    int solveMode=2;
1390    if (model->numberIterations())
1391      solveMode+=8;
1392    if (valuesPass)
1393      solveMode+= 4;
1394    coinFactorizationB_->setSolveMode(solveMode);
1395    while (status()<-98) {
1396     
1397      int i;
1398      int numberBasic=0;
1399      int numberRowBasic;
1400      // Move pivot variables across if they look good
1401      int * pivotTemp = model->rowArray(0)->getIndices();
1402      assert (!model->rowArray(0)->getNumElements());
1403      if (!matrix->rhsOffset(model)) {
1404        // Seems to prefer things in order so quickest
1405        // way is to go though like this
1406        for (i=0;i<numberRows;i++) {
1407          if (model->getRowStatus(i) == ClpSimplex::basic) 
1408            pivotTemp[numberBasic++]=i;
1409        }
1410        numberRowBasic=numberBasic;
1411        /* Put column basic variables into pivotVariable
1412           This is done by ClpMatrixBase to allow override for gub
1413        */
1414        matrix->generalExpanded(model,0,numberBasic);
1415      } else {
1416        // Long matrix - do a different way
1417        bool fullSearch=false;
1418        for (i=0;i<numberRows;i++) {
1419          int iPivot = pivotVariable[i];
1420          if (iPivot>=numberColumns) {
1421            pivotTemp[numberBasic++]=iPivot-numberColumns;
1422          }
1423        }
1424        numberRowBasic=numberBasic;
1425        for (i=0;i<numberRows;i++) {
1426          int iPivot = pivotVariable[i];
1427          if (iPivot<numberColumns) {
1428            if (iPivot>=0) {
1429              pivotTemp[numberBasic++]=iPivot;
1430            } else {
1431              // not full basis
1432              fullSearch=true;
1433              break;
1434            }
1435          }
1436        }
1437        if (fullSearch) {
1438          // do slow way
1439          numberBasic=0;
1440          for (i=0;i<numberRows;i++) {
1441            if (model->getRowStatus(i) == ClpSimplex::basic) 
1442              pivotTemp[numberBasic++]=i;
1443          }
1444          numberRowBasic=numberBasic;
1445          /* Put column basic variables into pivotVariable
1446             This is done by ClpMatrixBase to allow override for gub
1447          */
1448          matrix->generalExpanded(model,0,numberBasic);
1449        }
1450      }
1451      if (numberBasic>model->maximumBasic()) {
1452        // Take out some
1453        numberBasic=numberRowBasic;
1454        for (int i=0;i<numberColumns;i++) {
1455          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1456            if (numberBasic<numberRows)
1457              numberBasic++;
1458            else
1459              model->setColumnStatus(i,ClpSimplex::superBasic);
1460          }
1461        }
1462        numberBasic=numberRowBasic;
1463        matrix->generalExpanded(model,0,numberBasic);
1464      } else if (numberBasic<numberRows) {
1465        // add in some
1466        int needed = numberRows-numberBasic;
1467        // move up columns
1468        for (i=numberBasic-1;i>=numberRowBasic;i--)
1469          pivotTemp[i+needed]=pivotTemp[i];
1470        numberRowBasic=0;
1471        numberBasic=numberRows;
1472        for (i=0;i<numberRows;i++) {
1473          if (model->getRowStatus(i) == ClpSimplex::basic) {
1474            pivotTemp[numberRowBasic++]=i;
1475          } else if (needed) {
1476            needed--;
1477            model->setRowStatus(i,ClpSimplex::basic);
1478            pivotTemp[numberRowBasic++]=i;
1479          }
1480        }
1481      }
1482      CoinBigIndex numberElements=numberRowBasic;
1483     
1484      // compute how much in basis
1485      // can change for gub
1486      int numberColumnBasic = numberBasic-numberRowBasic;
1487     
1488      numberElements +=matrix->countBasis(model,
1489                                          pivotTemp+numberRowBasic, 
1490                                          numberRowBasic,
1491                                          numberColumnBasic);
1492      // Not needed for dense
1493      numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1494      int numberIterations=model->numberIterations();
1495      coinFactorizationB_->setUsefulInformation(&numberIterations,0);
1496      coinFactorizationB_->getAreas ( numberRows,
1497                                      numberRowBasic+numberColumnBasic, numberElements,
1498                                      2 * numberElements );
1499      // Fill in counts so we can skip part of preProcess
1500      // This is NOT needed for dense but would be needed for later versions
1501      CoinFactorizationDouble * elementU;
1502      int * indexRowU;
1503      CoinBigIndex * startColumnU;
1504      int * numberInRow;
1505      int * numberInColumn;
1506      elementU = coinFactorizationB_->elements();
1507      indexRowU = coinFactorizationB_->indices();
1508      startColumnU = coinFactorizationB_->starts();
1509#ifndef COIN_FAST_CODE
1510      double slackValue;
1511      slackValue = coinFactorizationB_->slackValue();
1512#else
1513#define slackValue -1.0
1514#endif
1515      numberInRow = coinFactorizationB_->numberInRow();
1516      numberInColumn = coinFactorizationB_->numberInColumn();
1517      CoinZeroN ( numberInRow, numberRows  );
1518      CoinZeroN ( numberInColumn, numberRows );
1519      for (i=0;i<numberRowBasic;i++) {
1520        int iRow = pivotTemp[i];
1521        // Change pivotTemp to correct sequence
1522        pivotTemp[i]=iRow+numberColumns;
1523        indexRowU[i]=iRow;
1524        startColumnU[i]=i;
1525        elementU[i]=slackValue;
1526        numberInRow[iRow]=1;
1527        numberInColumn[i]=1;
1528      }
1529      startColumnU[numberRowBasic]=numberRowBasic;
1530      // can change for gub so redo
1531      numberColumnBasic = numberBasic-numberRowBasic;
1532      matrix->fillBasis(model, 
1533                        pivotTemp+numberRowBasic, 
1534                        numberColumnBasic,
1535                        indexRowU, 
1536                        startColumnU+numberRowBasic,
1537                        numberInRow,
1538                        numberInColumn+numberRowBasic,
1539                        elementU);
1540      // recompute number basic
1541      numberBasic = numberRowBasic+numberColumnBasic;
1542      for (i=numberBasic;i<numberRows;i++)
1543        pivotTemp[i]=-1; // mark not there
1544      if (numberBasic) 
1545        numberElements = startColumnU[numberBasic-1]
1546          +numberInColumn[numberBasic-1];
1547      else
1548        numberElements=0;
1549      coinFactorizationB_->preProcess ( );
1550      coinFactorizationB_->factor (  );
1551      if (coinFactorizationB_->status() == -1 &&
1552          (coinFactorizationB_->solveMode()%3)!=0) {
1553        int solveMode=coinFactorizationB_->solveMode();
1554        solveMode -= solveMode%3; // so bottom will be 0
1555        coinFactorizationB_->setSolveMode(solveMode);
1556        coinFactorizationB_->setStatus(-99);
1557      }
1558      if (coinFactorizationB_->status() == -99)
1559        continue;
1560      // If we get here status is 0 or -1
1561      if (coinFactorizationB_->status() == 0&&numberBasic==numberRows) {
1562        coinFactorizationB_->postProcess(pivotTemp,pivotVariable);
1563      } else if (solveType==0||solveType==2/*||solveType==1*/) {
1564        // Change pivotTemp to be correct list
1565        anyChanged=true;
1566        coinFactorizationB_->makeNonSingular(pivotTemp,numberColumns);
1567        double * columnLower = model->lowerRegion();
1568        double * columnUpper = model->upperRegion();
1569        double * columnActivity = model->solutionRegion();
1570        double * rowLower = model->lowerRegion(0);
1571        double * rowUpper = model->upperRegion(0);
1572        double * rowActivity = model->solutionRegion(0);
1573        //redo basis - first take ALL out
1574        int iColumn;
1575        double largeValue = model->largeValue();
1576        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1577          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
1578            // take out
1579            if (!valuesPass) {
1580              double lower=columnLower[iColumn];
1581              double upper=columnUpper[iColumn];
1582              double value=columnActivity[iColumn];
1583              if (lower>-largeValue||upper<largeValue) {
1584                if (fabs(value-lower)<fabs(value-upper)) {
1585                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
1586                  columnActivity[iColumn]=lower;
1587                } else {
1588                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
1589                  columnActivity[iColumn]=upper;
1590                }
1591              } else {
1592                model->setColumnStatus(iColumn,ClpSimplex::isFree);
1593              }
1594            } else {
1595              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
1596            }
1597          }
1598        }
1599        int iRow;
1600        for (iRow=0;iRow<numberRows;iRow++) {
1601          if (model->getRowStatus(iRow)==ClpSimplex::basic) {
1602            // take out
1603            if (!valuesPass) {
1604              double lower=columnLower[iRow];
1605              double upper=columnUpper[iRow];
1606              double value=columnActivity[iRow];
1607              if (lower>-largeValue||upper<largeValue) {
1608                if (fabs(value-lower)<fabs(value-upper)) {
1609                  model->setRowStatus(iRow,ClpSimplex::atLowerBound);
1610                  columnActivity[iRow]=lower;
1611                } else {
1612                  model->setRowStatus(iRow,ClpSimplex::atUpperBound);
1613                  columnActivity[iRow]=upper;
1614                }
1615              } else {
1616                model->setRowStatus(iRow,ClpSimplex::isFree);
1617              }
1618            } else {
1619              model->setRowStatus(iRow,ClpSimplex::superBasic);
1620            }
1621          }
1622        }
1623        for (iRow=0;iRow<numberRows;iRow++) {
1624          int iSequence=pivotTemp[iRow];
1625          assert (iSequence>=0);
1626          // basic
1627          model->setColumnStatus(iSequence,ClpSimplex::basic);
1628        }
1629        // signal repeat
1630        coinFactorizationB_->setStatus(-99);
1631        // set fixed if they are
1632        for (iRow=0;iRow<numberRows;iRow++) {
1633          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
1634            if (rowLower[iRow]==rowUpper[iRow]) {
1635              rowActivity[iRow]=rowLower[iRow];
1636              model->setRowStatus(iRow,ClpSimplex::isFixed);
1637            }
1638          }
1639        }
1640        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1641          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
1642            if (columnLower[iColumn]==columnUpper[iColumn]) {
1643              columnActivity[iColumn]=columnLower[iColumn];
1644              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
1645            }
1646          }
1647        }
1648      } 
1649    }
1650#ifdef CLP_DEBUG
1651    // check basic
1652    CoinIndexedVector region1(2*numberRows);
1653    CoinIndexedVector region2B(2*numberRows);
1654    int iPivot;
1655    double * arrayB = region2B.denseVector();
1656    int i;
1657    for (iPivot=0;iPivot<numberRows;iPivot++) {
1658      int iSequence = pivotVariable[iPivot];
1659      model->unpack(&region2B,iSequence);
1660      coinFactorizationB_->updateColumn(&region1,&region2B);
1661      if (fabs(arrayB[iPivot]-1.0)<1.0e-4) {
1662        // OK?
1663        arrayB[iPivot]=0.0;
1664      } else {
1665        assert (fabs(arrayB[iPivot])<1.0e-4);
1666        for (i=0;i<numberRows;i++) {
1667          if (fabs(arrayB[i]-1.0)<1.0e-4)
1668            break;
1669        }
1670        assert (i<numberRows);
1671        printf("variable on row %d landed up on row %d\n",iPivot,i);
1672        arrayB[i]=0.0;
1673      }
1674      for (i=0;i<numberRows;i++)
1675        assert (fabs(arrayB[i])<1.0e-4);
1676      region2B.clear();
1677    }
1678#endif
1679#ifdef CLP_FACTORIZATION_INSTRUMENT
1680    factorization_instrument(2);
1681#endif
1682    if ( anyChanged&&model->algorithm()<0&&solveType>0) {
1683      double dummyCost;
1684      static_cast<ClpSimplexDual *> (model)->changeBounds(3,
1685                                                          NULL,dummyCost);
1686    }
1687    return coinFactorizationB_->status();
1688  }
1689  // If too many compressions increase area
1690  if (coinFactorizationA_->pivots()>1&&coinFactorizationA_->numberCompressions()*10 > coinFactorizationA_->pivots()+10) {
1691    coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor()* 1.1);
1692  }
1693  //int numberPivots=coinFactorizationA_->pivots();
1694#if 0
1695  if (model->algorithm()>0)
1696    numberSave=-1;
1697#endif
1698#ifndef SLIM_CLP
1699  if (!networkBasis_||doCheck) {
1700#endif
1701    coinFactorizationA_->setStatus(-99);
1702    int * pivotVariable = model->pivotVariable();
1703    int nTimesRound=0;
1704    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
1705    while (coinFactorizationA_->status()<-98) {
1706      nTimesRound++;
1707     
1708      int i;
1709      int numberBasic=0;
1710      int numberRowBasic;
1711      // Move pivot variables across if they look good
1712      int * pivotTemp = model->rowArray(0)->getIndices();
1713      assert (!model->rowArray(0)->getNumElements());
1714      if (!matrix->rhsOffset(model)) {
1715#if 0
1716        if (numberSave>0) {
1717          int nStill=0;
1718          int nAtBound=0;
1719          int nZeroDual=0;
1720          CoinIndexedVector * array = model->rowArray(3);
1721          CoinIndexedVector * objArray = model->columnArray(1);
1722          array->clear();
1723          objArray->clear();
1724          double * cost = model->costRegion();
1725          double tolerance = model->primalTolerance();
1726          double offset=0.0;
1727          for (i=0;i<numberRows;i++) {
1728            int iPivot = pivotVariable[i];
1729            if (iPivot<numberColumns&&isDense(iPivot)) {
1730              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
1731                nStill++;
1732                double value=model->solutionRegion()[iPivot];
1733                double dual = model->dualRowSolution()[i];
1734                double lower=model->lowerRegion()[iPivot];
1735                double upper=model->upperRegion()[iPivot];
1736                ClpSimplex::Status status;
1737                if (fabs(value-lower)<tolerance) {
1738                  status=ClpSimplex::atLowerBound;
1739                  nAtBound++;
1740                } else if (fabs(value-upper)<tolerance) {
1741                  nAtBound++;
1742                  status=ClpSimplex::atUpperBound;
1743                } else if (value>lower&&value<upper) {
1744                  status=ClpSimplex::superBasic;
1745                } else {
1746                  status=ClpSimplex::basic;
1747                }
1748                if (status!=ClpSimplex::basic) {
1749                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
1750                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
1751                    model->setRowStatus(i,ClpSimplex::basic);
1752                    pivotVariable[i]=i+numberColumns;
1753                    model->dualRowSolution()[i]=0.0;
1754                    model->djRegion(0)[i]=0.0;
1755                    array->add(i,dual);
1756                    offset += dual*model->solutionRegion(0)[i];
1757                  }
1758                }
1759                if (fabs(dual)<1.0e-5)
1760                  nZeroDual++;
1761              }
1762            }
1763          }
1764          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
1765                 numberSave,nStill,nAtBound,nZeroDual,offset);
1766          if (array->getNumElements()) {
1767            // modify costs
1768            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
1769                                               objArray);
1770            array->clear();
1771            int n=objArray->getNumElements();
1772            int * indices = objArray->getIndices();
1773            double * elements = objArray->denseVector();
1774            for (i=0;i<n;i++) {
1775              int iColumn = indices[i];
1776              cost[iColumn] -= elements[iColumn];
1777              elements[iColumn]=0.0;
1778            }
1779            objArray->setNumElements(0);
1780          }
1781        }
1782#endif
1783        // Seems to prefer things in order so quickest
1784        // way is to go though like this
1785        for (i=0;i<numberRows;i++) {
1786          if (model->getRowStatus(i) == ClpSimplex::basic) 
1787            pivotTemp[numberBasic++]=i;
1788        }
1789        numberRowBasic=numberBasic;
1790        /* Put column basic variables into pivotVariable
1791           This is done by ClpMatrixBase to allow override for gub
1792        */
1793        matrix->generalExpanded(model,0,numberBasic);
1794      } else {
1795        // Long matrix - do a different way
1796        bool fullSearch=false;
1797        for (i=0;i<numberRows;i++) {
1798          int iPivot = pivotVariable[i];
1799          if (iPivot>=numberColumns) {
1800            pivotTemp[numberBasic++]=iPivot-numberColumns;
1801          }
1802        }
1803        numberRowBasic=numberBasic;
1804        for (i=0;i<numberRows;i++) {
1805          int iPivot = pivotVariable[i];
1806          if (iPivot<numberColumns) {
1807            if (iPivot>=0) {
1808              pivotTemp[numberBasic++]=iPivot;
1809            } else {
1810              // not full basis
1811              fullSearch=true;
1812              break;
1813            }
1814          }
1815        }
1816        if (fullSearch) {
1817          // do slow way
1818          numberBasic=0;
1819          for (i=0;i<numberRows;i++) {
1820            if (model->getRowStatus(i) == ClpSimplex::basic) 
1821              pivotTemp[numberBasic++]=i;
1822          }
1823          numberRowBasic=numberBasic;
1824          /* Put column basic variables into pivotVariable
1825             This is done by ClpMatrixBase to allow override for gub
1826          */
1827          matrix->generalExpanded(model,0,numberBasic);
1828        }
1829      }
1830      if (numberBasic>model->maximumBasic()) {
1831#if 0 // ndef NDEBUG
1832        printf("%d basic - should only be %d\n",
1833               numberBasic,numberRows);
1834#endif
1835        // Take out some
1836        numberBasic=numberRowBasic;
1837        for (int i=0;i<numberColumns;i++) {
1838          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1839            if (numberBasic<numberRows)
1840              numberBasic++;
1841            else
1842              model->setColumnStatus(i,ClpSimplex::superBasic);
1843          }
1844        }
1845        numberBasic=numberRowBasic;
1846        matrix->generalExpanded(model,0,numberBasic);
1847      }
1848#ifndef SLIM_CLP
1849      // see if matrix a network
1850#ifndef NO_RTTI
1851      ClpNetworkMatrix* networkMatrix =
1852        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
1853#else
1854      ClpNetworkMatrix* networkMatrix = NULL;
1855      if (model->clpMatrix()->type()==11)
1856        networkMatrix = 
1857        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
1858#endif
1859      // If network - still allow ordinary factorization first time for laziness
1860      if (networkMatrix)
1861        coinFactorizationA_->setBiasLU(0); // All to U if network
1862      //int saveMaximumPivots = maximumPivots();
1863      delete networkBasis_;
1864      networkBasis_ = NULL;
1865      if (networkMatrix&&!doCheck)
1866        maximumPivots(1);
1867#endif
1868      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
1869      while (coinFactorizationA_->status()==-99) {
1870        // maybe for speed will be better to leave as many regions as possible
1871        coinFactorizationA_->gutsOfDestructor();
1872        coinFactorizationA_->gutsOfInitialize(2);
1873        CoinBigIndex numberElements=numberRowBasic;
1874
1875        // compute how much in basis
1876
1877        int i;
1878        // can change for gub
1879        int numberColumnBasic = numberBasic-numberRowBasic;
1880
1881        numberElements +=matrix->countBasis(model,
1882                                           pivotTemp+numberRowBasic, 
1883                                           numberRowBasic,
1884                                            numberColumnBasic);
1885        // and recompute as network side say different
1886        if (model->numberIterations())
1887          numberRowBasic = numberBasic - numberColumnBasic;
1888        numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1889        coinFactorizationA_->getAreas ( numberRows,
1890                   numberRowBasic+numberColumnBasic, numberElements,
1891                   2 * numberElements );
1892        //fill
1893        // Fill in counts so we can skip part of preProcess
1894        int * numberInRow = coinFactorizationA_->numberInRow();
1895        int * numberInColumn = coinFactorizationA_->numberInColumn();
1896        CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
1897        CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
1898        CoinFactorizationDouble * elementU = coinFactorizationA_->elementU();
1899        int * indexRowU = coinFactorizationA_->indexRowU();
1900        CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
1901#ifndef COIN_FAST_CODE
1902        double slackValue = coinFactorizationA_->slackValue();
1903#endif
1904        for (i=0;i<numberRowBasic;i++) {
1905          int iRow = pivotTemp[i];
1906          indexRowU[i]=iRow;
1907          startColumnU[i]=i;
1908          elementU[i]=slackValue;
1909          numberInRow[iRow]=1;
1910          numberInColumn[i]=1;
1911        }
1912        startColumnU[numberRowBasic]=numberRowBasic;
1913        // can change for gub so redo
1914        numberColumnBasic = numberBasic-numberRowBasic;
1915        matrix->fillBasis(model, 
1916                          pivotTemp+numberRowBasic, 
1917                          numberColumnBasic,
1918                          indexRowU, 
1919                          startColumnU+numberRowBasic,
1920                          numberInRow,
1921                          numberInColumn+numberRowBasic,
1922                          elementU);
1923#if 0
1924        {
1925          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
1926          for (int i=0;i<numberElements;i++)
1927            printf("row %d col %d value %g\n",indexRowU[i],indexColumnU_[i],
1928                   elementU[i]);
1929        }
1930#endif
1931        // recompute number basic
1932        numberBasic = numberRowBasic+numberColumnBasic;
1933        if (numberBasic) 
1934          numberElements = startColumnU[numberBasic-1]
1935            +numberInColumn[numberBasic-1];
1936        else
1937          numberElements=0;
1938        coinFactorizationA_->setNumberElementsU(numberElements);
1939        //saveFactorization("dump.d");
1940        if (coinFactorizationA_->biasLU()>=3||coinFactorizationA_->numberRows()!=coinFactorizationA_->numberColumns())
1941          coinFactorizationA_->preProcess ( 2 );
1942        else
1943          coinFactorizationA_->preProcess ( 3 ); // no row copy
1944        coinFactorizationA_->factor (  );
1945        if (coinFactorizationA_->status()==-99) {
1946          // get more memory
1947          coinFactorizationA_->areaFactor(2.0*coinFactorizationA_->areaFactor());
1948        } else if (coinFactorizationA_->status()==-1&&
1949                   (model->numberIterations()==0||nTimesRound>2)&&
1950                   coinFactorizationA_->denseThreshold()) {
1951          // Round again without dense
1952          coinFactorizationA_->setDenseThreshold(0);
1953          coinFactorizationA_->setStatus(-99);
1954        }
1955      }
1956      // If we get here status is 0 or -1
1957      if (coinFactorizationA_->status() == 0) {
1958        // We may need to tamper with order and redo - e.g. network with side
1959        int useNumberRows = numberRows;
1960        // **** we will also need to add test in dual steepest to do
1961        // as we do for network
1962        matrix->generalExpanded(model,12,useNumberRows);
1963        const int * permuteBack = coinFactorizationA_->permuteBack();
1964        const int * back = coinFactorizationA_->pivotColumnBack();
1965        //int * pivotTemp = pivotColumn_.array();
1966        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
1967        // Redo pivot order
1968        for (i=0;i<numberRowBasic;i++) {
1969          int k = pivotTemp[i];
1970          // so rowIsBasic[k] would be permuteBack[back[i]]
1971          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
1972        }
1973        for (;i<useNumberRows;i++) {
1974          int k = pivotTemp[i];
1975          // so rowIsBasic[k] would be permuteBack[back[i]]
1976          pivotVariable[permuteBack[back[i]]]=k;
1977        }
1978#if 0
1979        if (numberSave>=0) {
1980          numberSave=numberDense_;
1981          memset(saveList,0,((coinFactorizationA_->numberRows()+31)>>5)*sizeof(int));
1982          for (i=coinFactorizationA_->numberRows()-numberSave;i<coinFactorizationA_->numberRows();i++) {
1983            int k=pivotTemp[pivotColumn_.array()[i]];
1984            setDense(k);
1985          }
1986        }
1987#endif
1988        // Set up permutation vector
1989        // these arrays start off as copies of permute
1990        // (and we could use permute_ instead of pivotColumn (not back though))
1991        ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn()  );
1992        ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack()  );
1993#ifndef SLIM_CLP
1994        if (networkMatrix) {
1995          maximumPivots(CoinMax(2000,maximumPivots()));
1996          // redo arrays
1997          for (int iRow=0;iRow<4;iRow++) {
1998            int length =model->numberRows()+maximumPivots();
1999            if (iRow==3||model->objectiveAsObject()->type()>1)
2000              length += model->numberColumns();
2001            model->rowArray(iRow)->reserve(length);
2002          }
2003          // create network factorization
2004          if (doCheck)
2005            delete networkBasis_; // temp
2006          networkBasis_ = new ClpNetworkBasis(model,coinFactorizationA_->numberRows(),
2007                                              coinFactorizationA_->pivotRegion(),
2008                                              coinFactorizationA_->permuteBack(),
2009                                              coinFactorizationA_->startColumnU(),
2010                                              coinFactorizationA_->numberInColumn(),
2011                                              coinFactorizationA_->indexRowU(),
2012                                              coinFactorizationA_->elementU());
2013          // kill off arrays in ordinary factorization
2014          if (!doCheck) {
2015            coinFactorizationA_->gutsOfDestructor();
2016            // but make sure coinFactorizationA_->numberRows() set
2017            coinFactorizationA_->setNumberRows(model->numberRows());
2018            coinFactorizationA_->setStatus(0);
2019#if 0
2020            // but put back permute arrays so odd things will work
2021            int numberRows = model->numberRows();
2022            pivotColumnBack_ = new int [numberRows];
2023            permute_ = new int [numberRows];
2024            int i;
2025            for (i=0;i<numberRows;i++) {
2026              pivotColumnBack_[i]=i;
2027              permute_[i]=i;
2028            }
2029#endif
2030          }
2031        } else {
2032#endif
2033          // See if worth going sparse and when
2034          coinFactorizationA_->checkSparse();
2035#ifndef SLIM_CLP
2036        }
2037#endif
2038      } else if (coinFactorizationA_->status() == -1&&(solveType==0||solveType==2)) {
2039        // This needs redoing as it was merged coding - does not need array
2040        int numberTotal = numberRows+numberColumns;
2041        int * isBasic = new int [numberTotal];
2042        int * rowIsBasic = isBasic+numberColumns;
2043        int * columnIsBasic = isBasic;
2044        for (i=0;i<numberTotal;i++) 
2045          isBasic[i]=-1;
2046        for (i=0;i<numberRowBasic;i++) {
2047          int iRow = pivotTemp[i];
2048          rowIsBasic[iRow]=1;
2049        }
2050        for (;i<numberBasic;i++) {
2051          int iColumn = pivotTemp[i];
2052          columnIsBasic[iColumn]=1;
2053        }
2054        numberBasic=0;
2055        for (i=0;i<numberRows;i++) 
2056          pivotVariable[i]=-1;
2057        // mark as basic or non basic
2058        const int * pivotColumn = coinFactorizationA_->pivotColumn();
2059        for (i=0;i<numberRows;i++) {
2060          if (rowIsBasic[i]>=0) {
2061            if (pivotColumn[numberBasic]>=0) {
2062              rowIsBasic[i]=pivotColumn[numberBasic];
2063            } else {
2064              rowIsBasic[i]=-1;
2065              model->setRowStatus(i,ClpSimplex::superBasic);
2066            }
2067            numberBasic++;
2068          }
2069        }
2070        for (i=0;i<numberColumns;i++) {
2071          if (columnIsBasic[i]>=0) {
2072            if (pivotColumn[numberBasic]>=0) 
2073              columnIsBasic[i]=pivotColumn[numberBasic];
2074            else
2075              columnIsBasic[i]=-1;
2076            numberBasic++;
2077          }
2078        }
2079        // leave pivotVariable in useful form for cleaning basis
2080        int * pivotVariable = model->pivotVariable();
2081        for (i=0;i<numberRows;i++) {
2082          pivotVariable[i]=-1;
2083        }
2084
2085        for (i=0;i<numberRows;i++) {
2086          if (model->getRowStatus(i) == ClpSimplex::basic) {
2087            int iPivot = rowIsBasic[i];
2088            if (iPivot>=0) 
2089              pivotVariable[iPivot]=i+numberColumns;
2090          }
2091        }
2092        for (i=0;i<numberColumns;i++) {
2093          if (model->getColumnStatus(i) == ClpSimplex::basic) {
2094            int iPivot = columnIsBasic[i];
2095            if (iPivot>=0) 
2096              pivotVariable[iPivot]=i;
2097          }
2098        }
2099        delete [] isBasic;
2100        double * columnLower = model->lowerRegion();
2101        double * columnUpper = model->upperRegion();
2102        double * columnActivity = model->solutionRegion();
2103        double * rowLower = model->lowerRegion(0);
2104        double * rowUpper = model->upperRegion(0);
2105        double * rowActivity = model->solutionRegion(0);
2106        //redo basis - first take ALL columns out
2107        int iColumn;
2108        double largeValue = model->largeValue();
2109        for (iColumn=0;iColumn<numberColumns;iColumn++) {
2110          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
2111            // take out
2112            if (!valuesPass) {
2113              double lower=columnLower[iColumn];
2114              double upper=columnUpper[iColumn];
2115              double value=columnActivity[iColumn];
2116              if (lower>-largeValue||upper<largeValue) {
2117                if (fabs(value-lower)<fabs(value-upper)) {
2118                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
2119                  columnActivity[iColumn]=lower;
2120                } else {
2121                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
2122                  columnActivity[iColumn]=upper;
2123                }
2124              } else {
2125                model->setColumnStatus(iColumn,ClpSimplex::isFree);
2126              }
2127            } else {
2128              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
2129            }
2130          }
2131        }
2132        int iRow;
2133        for (iRow=0;iRow<numberRows;iRow++) {
2134          int iSequence=pivotVariable[iRow];
2135          if (iSequence>=0) {
2136            // basic
2137            if (iSequence>=numberColumns) {
2138              // slack in - leave
2139              //assert (iSequence-numberColumns==iRow);
2140            } else {
2141              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
2142              // put back structural
2143              model->setColumnStatus(iSequence,ClpSimplex::basic);
2144            }
2145          } else {
2146            // put in slack
2147            model->setRowStatus(iRow,ClpSimplex::basic);
2148          }
2149        }
2150        // Put back any key variables for gub
2151        int dummy;
2152        matrix->generalExpanded(model,1,dummy);
2153        // signal repeat
2154        coinFactorizationA_->setStatus(-99);
2155        // set fixed if they are
2156        for (iRow=0;iRow<numberRows;iRow++) {
2157          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
2158            if (rowLower[iRow]==rowUpper[iRow]) {
2159              rowActivity[iRow]=rowLower[iRow];
2160              model->setRowStatus(iRow,ClpSimplex::isFixed);
2161            }
2162          }
2163        }
2164        for (iColumn=0;iColumn<numberColumns;iColumn++) {
2165          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
2166            if (columnLower[iColumn]==columnUpper[iColumn]) {
2167              columnActivity[iColumn]=columnLower[iColumn];
2168              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
2169            }
2170          }
2171        }
2172      } 
2173    }
2174#ifndef SLIM_CLP
2175  } else {
2176    // network - fake factorization - do nothing
2177    coinFactorizationA_->setStatus(0);
2178    coinFactorizationA_->setPivots(0);
2179  }
2180#endif
2181#ifndef SLIM_CLP
2182  if (!coinFactorizationA_->status()) {
2183    // take out part if quadratic
2184    if (model->algorithm()==2) {
2185      ClpObjective * obj = model->objectiveAsObject();
2186#ifndef NDEBUG
2187      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2188      assert (quadraticObj);
2189#else
2190      ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
2191#endif
2192      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2193      int numberXColumns = quadratic->getNumCols();
2194      assert (numberXColumns<numberColumns);
2195      int base = numberColumns-numberXColumns;
2196      int * which = new int [numberXColumns];
2197      int * pivotVariable = model->pivotVariable();
2198      int * permute = pivotColumn();
2199      int i;
2200      int n=0;
2201      for (i=0;i<numberRows;i++) {
2202        int iSj = pivotVariable[i]-base;
2203        if (iSj>=0&&iSj<numberXColumns) 
2204          which[n++]=permute[i];
2205      }
2206      if (n)
2207        coinFactorizationA_->emptyRows(n,which);
2208      delete [] which;
2209    }
2210  }
2211#endif
2212#ifdef CLP_FACTORIZATION_INSTRUMENT
2213  factorization_instrument(2);
2214#endif
2215  return coinFactorizationA_->status();
2216}
2217/* Replaces one Column in basis,
2218   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
2219   If checkBeforeModifying is true will do all accuracy checks
2220   before modifying factorization.  Whether to set this depends on
2221   speed considerations.  You could just do this on first iteration
2222   after factorization and thereafter re-factorize
2223   partial update already in U */
2224int 
2225ClpFactorization::replaceColumn ( const ClpSimplex * model, 
2226                                  CoinIndexedVector * regionSparse,
2227                                  CoinIndexedVector * tableauColumn,
2228                                  int pivotRow,
2229                                  double pivotCheck ,
2230                                  bool checkBeforeModifying)
2231{
2232#ifndef SLIM_CLP
2233  if (!networkBasis_) {
2234#endif
2235#ifdef CLP_FACTORIZATION_INSTRUMENT
2236    factorization_instrument(-1);
2237#endif
2238    int returnCode;
2239    // see if FT
2240    if (!coinFactorizationA_||coinFactorizationA_->forrestTomlin()) {
2241      if (coinFactorizationA_) {
2242        returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2243                                                        pivotRow,
2244                                                        pivotCheck,
2245                                                        checkBeforeModifying);
2246      } else {
2247        bool tab = coinFactorizationB_->wantsTableauColumn();
2248        int numberIterations=model->numberIterations();
2249        coinFactorizationB_->setUsefulInformation(&numberIterations,1);
2250        returnCode= 
2251          coinFactorizationB_->replaceColumn(tab?tableauColumn:regionSparse,
2252                                             pivotRow,
2253                                             pivotCheck,
2254                                             checkBeforeModifying);
2255#ifdef CLP_DEBUG
2256        // check basic
2257        int numberRows = coinFactorizationB_->numberRows();
2258        CoinIndexedVector region1(2*numberRows);
2259        CoinIndexedVector region2A(2*numberRows);
2260        CoinIndexedVector region2B(2*numberRows);
2261        int iPivot;
2262        double * arrayB = region2B.denseVector();
2263        int * pivotVariable = model->pivotVariable();
2264        int i;
2265        for (iPivot=0;iPivot<numberRows;iPivot++) {
2266          int iSequence = pivotVariable[iPivot];
2267          if (iPivot==pivotRow)
2268            iSequence = model->sequenceIn();
2269          model->unpack(&region2B,iSequence);
2270          coinFactorizationB_->updateColumn(&region1,&region2B);
2271          assert (fabs(arrayB[iPivot]-1.0)<1.0e-4);
2272          arrayB[iPivot]=0.0;
2273          for (i=0;i<numberRows;i++)
2274            assert (fabs(arrayB[i])<1.0e-4);
2275          region2B.clear();
2276        }
2277#endif
2278      }
2279    } else {
2280      returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn,
2281                                              pivotRow,pivotCheck); // Note array
2282    }
2283#ifdef CLP_FACTORIZATION_INSTRUMENT
2284      factorization_instrument(3);
2285#endif
2286      return returnCode;
2287
2288#ifndef SLIM_CLP
2289  } else {
2290    if (doCheck) {
2291      int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2292                                                pivotRow,
2293                                                pivotCheck,
2294                                                checkBeforeModifying);
2295      networkBasis_->replaceColumn(regionSparse,
2296                                   pivotRow);
2297      return returnCode;
2298    } else {
2299      // increase number of pivots
2300      coinFactorizationA_->setPivots(coinFactorizationA_->pivots()+1);
2301      return networkBasis_->replaceColumn(regionSparse,
2302                                   pivotRow);
2303    }
2304  }
2305#endif
2306}
2307
2308/* Updates one column (FTRAN) from region2
2309   number returned is negative if no room
2310   region1 starts as zero and is zero at end */
2311int 
2312ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
2313                                   CoinIndexedVector * regionSparse2)
2314{
2315#ifdef CLP_DEBUG
2316  regionSparse->checkClear();
2317#endif
2318  if (!numberRows())
2319      return 0;
2320#ifndef SLIM_CLP
2321  if (!networkBasis_) {
2322#endif
2323#ifdef CLP_FACTORIZATION_INSTRUMENT
2324    factorization_instrument(-1);
2325#endif
2326    int returnCode;
2327    if (coinFactorizationA_) {
2328      coinFactorizationA_->setCollectStatistics(true);
2329      returnCode= coinFactorizationA_->updateColumnFT(regionSparse,
2330                                                           regionSparse2);
2331      coinFactorizationA_->setCollectStatistics(false);
2332    } else {
2333      returnCode= coinFactorizationB_->updateColumnFT(regionSparse,
2334                                                         regionSparse2);
2335    }
2336#ifdef CLP_FACTORIZATION_INSTRUMENT
2337    factorization_instrument(4);
2338#endif
2339    return returnCode;
2340#ifndef SLIM_CLP
2341  } else {
2342#ifdef CHECK_NETWORK
2343    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2344    double * check = new double[coinFactorizationA_->numberRows()];
2345    int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2346                                                       regionSparse2);
2347    networkBasis_->updateColumn(regionSparse,save,-1);
2348    int i;
2349    double * array = regionSparse2->denseVector();
2350    int * indices = regionSparse2->getIndices();
2351    int n=regionSparse2->getNumElements();
2352    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2353    double * array2 = save->denseVector();
2354    int * indices2 = save->getIndices();
2355    int n2=save->getNumElements();
2356    assert (n==n2);
2357    if (save->packedMode()) {
2358      for (i=0;i<n;i++) {
2359        check[indices[i]]=array[i];
2360      }
2361      for (i=0;i<n;i++) {
2362        double value2 = array2[i];
2363        assert (check[indices2[i]]==value2);
2364      }
2365    } else {
2366      int numberRows = coinFactorizationA_->numberRows();
2367      for (i=0;i<numberRows;i++) {
2368        double value1 = array[i];
2369        double value2 = array2[i];
2370        assert (value1==value2);
2371      }
2372    }
2373    delete save;
2374    delete [] check;
2375    return returnCode;
2376#else
2377    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
2378    return 1;
2379#endif
2380  }
2381#endif
2382}
2383/* Updates one column (FTRAN) from region2
2384   number returned is negative if no room
2385   region1 starts as zero and is zero at end */
2386int 
2387ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
2388                                 CoinIndexedVector * regionSparse2,
2389                                 bool noPermute) const
2390{
2391#ifdef CLP_DEBUG
2392  if (!noPermute)
2393    regionSparse->checkClear();
2394#endif
2395  if (!numberRows())
2396    return 0;
2397#ifndef SLIM_CLP
2398  if (!networkBasis_) {
2399#endif
2400#ifdef CLP_FACTORIZATION_INSTRUMENT
2401    factorization_instrument(-1);
2402#endif
2403    int returnCode;
2404    if (coinFactorizationA_) {
2405      coinFactorizationA_->setCollectStatistics(true);
2406      returnCode= coinFactorizationA_->updateColumn(regionSparse,
2407                                                         regionSparse2,
2408                                                         noPermute);
2409      coinFactorizationA_->setCollectStatistics(false);
2410    } else {
2411      returnCode= coinFactorizationB_->updateColumn(regionSparse,
2412                                                         regionSparse2,
2413                                                         noPermute);
2414    }
2415#ifdef CLP_FACTORIZATION_INSTRUMENT
2416    factorization_instrument(5);
2417#endif
2418    //#define PRINT_VECTOR
2419#ifdef PRINT_VECTOR
2420    printf("Update\n");
2421    regionSparse2->print();
2422#endif
2423    return returnCode;
2424#ifndef SLIM_CLP
2425  } else {
2426#ifdef CHECK_NETWORK
2427    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2428    double * check = new double[coinFactorizationA_->numberRows()];
2429    int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2430                                                     regionSparse2,
2431                                                     noPermute);
2432    networkBasis_->updateColumn(regionSparse,save,-1);
2433    int i;
2434    double * array = regionSparse2->denseVector();
2435    int * indices = regionSparse2->getIndices();
2436    int n=regionSparse2->getNumElements();
2437    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2438    double * array2 = save->denseVector();
2439    int * indices2 = save->getIndices();
2440    int n2=save->getNumElements();
2441    assert (n==n2);
2442    if (save->packedMode()) {
2443      for (i=0;i<n;i++) {
2444        check[indices[i]]=array[i];
2445      }
2446      for (i=0;i<n;i++) {
2447        double value2 = array2[i];
2448        assert (check[indices2[i]]==value2);
2449      }
2450    } else {
2451      int numberRows = coinFactorizationA_->numberRows();
2452      for (i=0;i<numberRows;i++) {
2453        double value1 = array[i];
2454        double value2 = array2[i];
2455        assert (value1==value2);
2456      }
2457    }
2458    delete save;
2459    delete [] check;
2460    return returnCode;
2461#else
2462    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
2463    return 1;
2464#endif
2465  }
2466#endif
2467}
2468/* Updates one column (FTRAN) from region2
2469   Tries to do FT update
2470   number returned is negative if no room.
2471   Also updates region3
2472   region1 starts as zero and is zero at end */
2473int 
2474ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
2475                                       CoinIndexedVector * regionSparse2,
2476                                       CoinIndexedVector * regionSparse3,
2477                                       bool noPermuteRegion3)
2478{
2479#ifdef CLP_DEBUG
2480  regionSparse1->checkClear();
2481#endif
2482  if (!numberRows())
2483    return 0;
2484  int returnCode=0;
2485#ifndef SLIM_CLP
2486  if (!networkBasis_) {
2487#endif
2488#ifdef CLP_FACTORIZATION_INSTRUMENT
2489    factorization_instrument(-1);
2490#endif
2491    if (coinFactorizationA_) {
2492      coinFactorizationA_->setCollectStatistics(true);
2493      if (coinFactorizationA_->spaceForForrestTomlin()) {
2494        assert (regionSparse2->packedMode());
2495        assert (!regionSparse3->packedMode());
2496        returnCode= coinFactorizationA_->updateTwoColumnsFT(regionSparse1,
2497                                                            regionSparse2,
2498                                                            regionSparse3,
2499                                                            noPermuteRegion3);
2500      } else {
2501        returnCode= coinFactorizationA_->updateColumnFT(regionSparse1,
2502                                                        regionSparse2);
2503        coinFactorizationA_->updateColumn(regionSparse1,
2504                                          regionSparse3,
2505                                          noPermuteRegion3);
2506      }
2507      coinFactorizationA_->setCollectStatistics(false);
2508    } else {
2509#if 0
2510      CoinSimpFactorization * fact =
2511        dynamic_cast< CoinSimpFactorization*>(coinFactorizationB_);
2512      if (!fact) {
2513        returnCode= coinFactorizationB_->updateColumnFT(regionSparse1,
2514                                                        regionSparse2);
2515        coinFactorizationB_->updateColumn(regionSparse1,
2516                                          regionSparse3,
2517                                          noPermuteRegion3);
2518      } else {
2519        returnCode= fact->updateTwoColumnsFT(regionSparse1,
2520                                             regionSparse2,
2521                                             regionSparse3,
2522                                             noPermuteRegion3);
2523      }
2524#else
2525      returnCode= coinFactorizationB_->updateTwoColumnsFT(regionSparse1,
2526                                                          regionSparse2,
2527                                                          regionSparse3,
2528                                                          noPermuteRegion3);
2529#endif
2530    }
2531#ifdef CLP_FACTORIZATION_INSTRUMENT
2532    factorization_instrument(9);
2533#endif
2534#ifdef PRINT_VECTOR
2535    printf("UpdateTwoFT\n");
2536    regionSparse2->print();
2537    regionSparse3->print();
2538#endif
2539    return returnCode;
2540#ifndef SLIM_CLP
2541  } else {
2542    returnCode=updateColumnFT(regionSparse1,regionSparse2);
2543    updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
2544  }
2545#endif
2546  return returnCode;
2547}
2548/* Updates one column (FTRAN) from region2
2549   number returned is negative if no room
2550   region1 starts as zero and is zero at end */
2551int 
2552ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
2553                                 CoinIndexedVector * regionSparse2,
2554                                 bool noPermute) const
2555{
2556  if (!noPermute)
2557    regionSparse->checkClear();
2558  if (!coinFactorizationA_->numberRows())
2559    return 0;
2560  coinFactorizationA_->setCollectStatistics(false);
2561  int returnCode= coinFactorizationA_->updateColumn(regionSparse,
2562                                                   regionSparse2,
2563                                                   noPermute);
2564  return returnCode;
2565}
2566/* Updates one column (BTRAN) from region2
2567   region1 starts as zero and is zero at end */
2568int 
2569ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
2570                          CoinIndexedVector * regionSparse2) const
2571{
2572  if (!numberRows())
2573    return 0;
2574#ifndef SLIM_CLP
2575  if (!networkBasis_) {
2576#endif
2577#ifdef CLP_FACTORIZATION_INSTRUMENT
2578    factorization_instrument(-1);
2579#endif
2580    int returnCode;
2581
2582    if (coinFactorizationA_) {
2583      coinFactorizationA_->setCollectStatistics(true);
2584      returnCode =  coinFactorizationA_->updateColumnTranspose(regionSparse,
2585                                                    regionSparse2);
2586      coinFactorizationA_->setCollectStatistics(false);
2587    } else {
2588      returnCode= coinFactorizationB_->updateColumnTranspose(regionSparse,
2589                                                         regionSparse2);
2590    }
2591#ifdef CLP_FACTORIZATION_INSTRUMENT
2592    factorization_instrument(6);
2593#endif
2594#ifdef PRINT_VECTOR
2595    printf("UpdateTranspose\n");
2596    regionSparse2->print();
2597#endif
2598    return returnCode;
2599#ifndef SLIM_CLP
2600  } else {
2601#ifdef CHECK_NETWORK
2602    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2603    double * check = new double[coinFactorizationA_->numberRows()];
2604    int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
2605                                                              regionSparse2);
2606    networkBasis_->updateColumnTranspose(regionSparse,save);
2607    int i;
2608    double * array = regionSparse2->denseVector();
2609    int * indices = regionSparse2->getIndices();
2610    int n=regionSparse2->getNumElements();
2611    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2612    double * array2 = save->denseVector();
2613    int * indices2 = save->getIndices();
2614    int n2=save->getNumElements();
2615    assert (n==n2);
2616    if (save->packedMode()) {
2617      for (i=0;i<n;i++) {
2618        check[indices[i]]=array[i];
2619      }
2620      for (i=0;i<n;i++) {
2621        double value2 = array2[i];
2622        assert (check[indices2[i]]==value2);
2623      }
2624    } else {
2625      int numberRows = coinFactorizationA_->numberRows();
2626      for (i=0;i<numberRows;i++) {
2627        double value1 = array[i];
2628        double value2 = array2[i];
2629        assert (value1==value2);
2630      }
2631    }
2632    delete save;
2633    delete [] check;
2634    return returnCode;
2635#else
2636    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
2637#endif
2638  }
2639#endif
2640}
2641/* makes a row copy of L for speed and to allow very sparse problems */
2642void 
2643ClpFactorization::goSparse()
2644{
2645#ifndef SLIM_CLP
2646  if (!networkBasis_) {
2647#endif
2648    if (coinFactorizationA_) {
2649#ifdef CLP_FACTORIZATION_INSTRUMENT
2650      factorization_instrument(-1);
2651#endif
2652      coinFactorizationA_->goSparse();
2653#ifdef CLP_FACTORIZATION_INSTRUMENT
2654      factorization_instrument(7);
2655#endif
2656    }
2657  }
2658}
2659// Cleans up i.e. gets rid of network basis
2660void 
2661ClpFactorization::cleanUp()
2662{
2663#ifndef SLIM_CLP
2664  delete networkBasis_;
2665  networkBasis_=NULL;
2666#endif
2667  if (coinFactorizationA_)
2668  coinFactorizationA_->resetStatistics();
2669}
2670/// Says whether to redo pivot order
2671bool 
2672ClpFactorization::needToReorder() const
2673{
2674#ifdef CHECK_NETWORK
2675  return true;
2676#endif
2677#ifndef SLIM_CLP
2678  if (!networkBasis_)
2679#endif
2680    return true;
2681#ifndef SLIM_CLP
2682  else
2683    return false;
2684#endif
2685}
2686// Get weighted row list
2687void
2688ClpFactorization::getWeights(int * weights) const
2689{
2690#ifdef CLP_FACTORIZATION_INSTRUMENT
2691  factorization_instrument(-1);
2692#endif
2693#ifndef SLIM_CLP
2694  if (networkBasis_) {
2695    // Network - just unit
2696    int numberRows = coinFactorizationA_->numberRows();
2697    for (int i=0;i<numberRows;i++) 
2698      weights[i]=1;
2699    return;
2700  }
2701#endif
2702  int * numberInRow = coinFactorizationA_->numberInRow();
2703  int * numberInColumn = coinFactorizationA_->numberInColumn();
2704  int * permuteBack = coinFactorizationA_->pivotColumnBack();
2705  int * indexRowU = coinFactorizationA_->indexRowU();
2706  const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
2707  const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
2708  int numberRows = coinFactorizationA_->numberRows();
2709  if (!startRowL||!coinFactorizationA_->numberInRow()) {
2710    int * temp = new int[numberRows];
2711    memset(temp,0,numberRows*sizeof(int));
2712    int i;
2713    for (i=0;i<numberRows;i++) {
2714      // one for pivot
2715      temp[i]++;
2716      CoinBigIndex j;
2717      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
2718        int iRow=indexRowU[j];
2719        temp[iRow]++;
2720      }
2721    }
2722    CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
2723    int * indexRowL = coinFactorizationA_->indexRowL();
2724    int numberL = coinFactorizationA_->numberL();
2725    CoinBigIndex baseL = coinFactorizationA_->baseL();
2726    for (i=baseL;i<baseL+numberL;i++) {
2727      CoinBigIndex j;
2728      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
2729        int iRow = indexRowL[j];
2730        temp[iRow]++;
2731      }
2732    }
2733    for (i=0;i<numberRows;i++) {
2734      int number = temp[i];
2735      int iPermute = permuteBack[i];
2736      weights[iPermute]=number;
2737    }
2738    delete [] temp;
2739  } else {
2740    int i;
2741    for (i=0;i<numberRows;i++) {
2742      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
2743      //number = startRowL[i+1]-startRowL[i]+1;
2744      //number = numberInRow[i]+1;
2745      int iPermute = permuteBack[i];
2746      weights[iPermute]=number;
2747    }
2748  }
2749#ifdef CLP_FACTORIZATION_INSTRUMENT
2750  factorization_instrument(8);
2751#endif
2752}
2753// Set tolerances to safer of existing and given
2754void 
2755ClpFactorization::saferTolerances (  double zeroValue, 
2756                                    double pivotValue)
2757{
2758  // better to have small tolerance even if slower
2759  zeroTolerance(CoinMin(zeroTolerance(),zeroValue));
2760  // better to have large tolerance even if slower
2761  double newValue;
2762  if (pivotValue>0.0)
2763    newValue = pivotValue;
2764  else 
2765    newValue = -pivotTolerance()*pivotValue;
2766  pivotTolerance(CoinMin(CoinMax(pivotTolerance(),newValue),0.999));
2767}
2768// Sets factorization
2769void 
2770ClpFactorization::setFactorization(ClpFactorization & rhs)
2771{
2772  ClpFactorization::operator=(rhs);
2773}
2774#endif
Note: See TracBrowser for help on using the repository browser.