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

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

mistake on osl

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 85.3 KB
Line 
1/* $Id: ClpFactorization.cpp 1395 2009-07-09 14:50:07Z 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                                  double acceptablePivot)
727{
728  int returnCode;
729#ifndef SLIM_CLP
730  if (!networkBasis_) {
731#endif
732#ifdef CLP_FACTORIZATION_INSTRUMENT
733    factorization_instrument(-1);
734#endif
735    // see if FT
736    if (doForrestTomlin_) {
737      returnCode= CoinFactorization::replaceColumn(regionSparse,
738                                              pivotRow,
739                                              pivotCheck,
740                                                   checkBeforeModifying,
741                                                   acceptablePivot);
742    } else {
743      returnCode= CoinFactorization::replaceColumnPFI(tableauColumn,
744                                              pivotRow,pivotCheck); // Note array
745    }
746#ifdef CLP_FACTORIZATION_INSTRUMENT
747    factorization_instrument(3);
748#endif
749
750#ifndef SLIM_CLP
751  } else {
752    if (doCheck) {
753      returnCode = CoinFactorization::replaceColumn(regionSparse,
754                                                        pivotRow,
755                                                        pivotCheck,
756                                                    checkBeforeModifying,
757                                                    acceptablePivot);
758      networkBasis_->replaceColumn(regionSparse,
759                                   pivotRow);
760    } else {
761      // increase number of pivots
762      numberPivots_++;
763      returnCode = networkBasis_->replaceColumn(regionSparse,
764                                   pivotRow);
765    }
766  }
767#endif
768  return returnCode;
769}
770
771/* Updates one column (FTRAN) from region2
772   number returned is negative if no room
773   region1 starts as zero and is zero at end */
774int 
775ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
776                                   CoinIndexedVector * regionSparse2)
777{
778#ifdef CLP_DEBUG
779  regionSparse->checkClear();
780#endif
781  if (!numberRows_)
782    return 0;
783#ifndef SLIM_CLP
784  if (!networkBasis_) {
785#endif
786#ifdef CLP_FACTORIZATION_INSTRUMENT
787    factorization_instrument(-1);
788#endif
789    collectStatistics_ = true;
790    int returnCode= CoinFactorization::updateColumnFT(regionSparse,
791                                                       regionSparse2);
792    collectStatistics_ = false;
793#ifdef CLP_FACTORIZATION_INSTRUMENT
794    factorization_instrument(4);
795#endif
796    return returnCode;
797#ifndef SLIM_CLP
798  } else {
799#ifdef CHECK_NETWORK
800    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
801    double * check = new double[numberRows_];
802    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
803                                                       regionSparse2);
804    networkBasis_->updateColumn(regionSparse,save,-1);
805    int i;
806    double * array = regionSparse2->denseVector();
807    int * indices = regionSparse2->getIndices();
808    int n=regionSparse2->getNumElements();
809    memset(check,0,numberRows_*sizeof(double));
810    double * array2 = save->denseVector();
811    int * indices2 = save->getIndices();
812    int n2=save->getNumElements();
813    assert (n==n2);
814    if (save->packedMode()) {
815      for (i=0;i<n;i++) {
816        check[indices[i]]=array[i];
817      }
818      for (i=0;i<n;i++) {
819        double value2 = array2[i];
820        assert (check[indices2[i]]==value2);
821      }
822    } else {
823      for (i=0;i<numberRows_;i++) {
824        double value1 = array[i];
825        double value2 = array2[i];
826        assert (value1==value2);
827      }
828    }
829    delete save;
830    delete [] check;
831    return returnCode;
832#else
833    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
834    return 1;
835#endif
836  }
837#endif
838}
839/* Updates one column (FTRAN) from region2
840   number returned is negative if no room
841   region1 starts as zero and is zero at end */
842int 
843ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
844                                 CoinIndexedVector * regionSparse2,
845                                 bool noPermute) const
846{
847#ifdef CLP_DEBUG
848  if (!noPermute)
849    regionSparse->checkClear();
850#endif
851  if (!numberRows_)
852    return 0;
853#ifndef SLIM_CLP
854  if (!networkBasis_) {
855#endif
856#ifdef CLP_FACTORIZATION_INSTRUMENT
857    factorization_instrument(-1);
858#endif
859    collectStatistics_ = true;
860    int returnCode= CoinFactorization::updateColumn(regionSparse,
861                                                     regionSparse2,
862                                                     noPermute);
863    collectStatistics_ = false;
864#ifdef CLP_FACTORIZATION_INSTRUMENT
865    factorization_instrument(5);
866#endif
867    return returnCode;
868#ifndef SLIM_CLP
869  } else {
870#ifdef CHECK_NETWORK
871    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
872    double * check = new double[numberRows_];
873    int returnCode = CoinFactorization::updateColumn(regionSparse,
874                                                     regionSparse2,
875                                                     noPermute);
876    networkBasis_->updateColumn(regionSparse,save,-1);
877    int i;
878    double * array = regionSparse2->denseVector();
879    int * indices = regionSparse2->getIndices();
880    int n=regionSparse2->getNumElements();
881    memset(check,0,numberRows_*sizeof(double));
882    double * array2 = save->denseVector();
883    int * indices2 = save->getIndices();
884    int n2=save->getNumElements();
885    assert (n==n2);
886    if (save->packedMode()) {
887      for (i=0;i<n;i++) {
888        check[indices[i]]=array[i];
889      }
890      for (i=0;i<n;i++) {
891        double value2 = array2[i];
892        assert (check[indices2[i]]==value2);
893      }
894    } else {
895      for (i=0;i<numberRows_;i++) {
896        double value1 = array[i];
897        double value2 = array2[i];
898        assert (value1==value2);
899      }
900    }
901    delete save;
902    delete [] check;
903    return returnCode;
904#else
905    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
906    return 1;
907#endif
908  }
909#endif
910}
911/* Updates one column (FTRAN) from region2
912   Tries to do FT update
913   number returned is negative if no room.
914   Also updates region3
915   region1 starts as zero and is zero at end */
916int 
917ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
918                                       CoinIndexedVector * regionSparse2,
919                                       CoinIndexedVector * regionSparse3,
920                                       bool noPermuteRegion3) 
921{
922  int returnCode=updateColumnFT(regionSparse1,regionSparse2);
923  updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
924  return returnCode;
925}
926/* Updates one column (FTRAN) from region2
927   number returned is negative if no room
928   region1 starts as zero and is zero at end */
929int 
930ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
931                                 CoinIndexedVector * regionSparse2,
932                                 bool noPermute) const
933{
934  if (!noPermute)
935    regionSparse->checkClear();
936  if (!numberRows_)
937    return 0;
938  collectStatistics_ = false;
939  int returnCode= CoinFactorization::updateColumn(regionSparse,
940                                                   regionSparse2,
941                                                   noPermute);
942  return returnCode;
943}
944/* Updates one column (BTRAN) from region2
945   region1 starts as zero and is zero at end */
946int 
947ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
948                          CoinIndexedVector * regionSparse2) const
949{
950  if (!numberRows_)
951    return 0;
952#ifndef SLIM_CLP
953  if (!networkBasis_) {
954#endif
955#ifdef CLP_FACTORIZATION_INSTRUMENT
956    factorization_instrument(-1);
957#endif
958    collectStatistics_ = true;
959    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
960                                                    regionSparse2);
961    collectStatistics_ = false;
962#ifdef CLP_FACTORIZATION_INSTRUMENT
963    factorization_instrument(6);
964#endif
965    return returnCode;
966#ifndef SLIM_CLP
967  } else {
968#ifdef CHECK_NETWORK
969    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
970    double * check = new double[numberRows_];
971    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
972                                                              regionSparse2);
973    networkBasis_->updateColumnTranspose(regionSparse,save);
974    int i;
975    double * array = regionSparse2->denseVector();
976    int * indices = regionSparse2->getIndices();
977    int n=regionSparse2->getNumElements();
978    memset(check,0,numberRows_*sizeof(double));
979    double * array2 = save->denseVector();
980    int * indices2 = save->getIndices();
981    int n2=save->getNumElements();
982    assert (n==n2);
983    if (save->packedMode()) {
984      for (i=0;i<n;i++) {
985        check[indices[i]]=array[i];
986      }
987      for (i=0;i<n;i++) {
988        double value2 = array2[i];
989        assert (check[indices2[i]]==value2);
990      }
991    } else {
992      for (i=0;i<numberRows_;i++) {
993        double value1 = array[i];
994        double value2 = array2[i];
995        assert (value1==value2);
996      }
997    }
998    delete save;
999    delete [] check;
1000    return returnCode;
1001#else
1002    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
1003#endif
1004  }
1005#endif
1006}
1007/* makes a row copy of L for speed and to allow very sparse problems */
1008void 
1009ClpFactorization::goSparse()
1010{
1011#ifdef CLP_FACTORIZATION_INSTRUMENT
1012  factorization_instrument(-1);
1013#endif
1014#ifndef SLIM_CLP
1015  if (!networkBasis_) 
1016#endif
1017    CoinFactorization::goSparse();
1018#ifdef CLP_FACTORIZATION_INSTRUMENT
1019  factorization_instrument(7);
1020#endif
1021}
1022// Cleans up i.e. gets rid of network basis
1023void 
1024ClpFactorization::cleanUp()
1025{
1026#ifndef SLIM_CLP
1027  delete networkBasis_;
1028  networkBasis_=NULL;
1029#endif
1030  resetStatistics();
1031}
1032/// Says whether to redo pivot order
1033bool 
1034ClpFactorization::needToReorder() const
1035{
1036#ifdef CHECK_NETWORK
1037  return true;
1038#endif
1039#ifndef SLIM_CLP
1040  if (!networkBasis_)
1041#endif
1042    return true;
1043#ifndef SLIM_CLP
1044  else
1045    return false;
1046#endif
1047}
1048// Get weighted row list
1049void
1050ClpFactorization::getWeights(int * weights) const
1051{
1052#ifdef CLP_FACTORIZATION_INSTRUMENT
1053  factorization_instrument(-1);
1054#endif
1055#ifndef SLIM_CLP
1056  if (networkBasis_) {
1057    // Network - just unit
1058    for (int i=0;i<numberRows_;i++) 
1059      weights[i]=1;
1060    return;
1061  }
1062#endif
1063  int * numberInRow = numberInRow_.array();
1064  int * numberInColumn = numberInColumn_.array();
1065  int * permuteBack = pivotColumnBack_.array();
1066  int * indexRowU = indexRowU_.array();
1067  const CoinBigIndex * startColumnU = startColumnU_.array();
1068  const CoinBigIndex * startRowL = startRowL_.array();
1069  if (!startRowL||!numberInRow_.array()) {
1070    int * temp = new int[numberRows_];
1071    memset(temp,0,numberRows_*sizeof(int));
1072    int i;
1073    for (i=0;i<numberRows_;i++) {
1074      // one for pivot
1075      temp[i]++;
1076      CoinBigIndex j;
1077      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
1078        int iRow=indexRowU[j];
1079        temp[iRow]++;
1080      }
1081    }
1082    CoinBigIndex * startColumnL = startColumnL_.array();
1083    int * indexRowL = indexRowL_.array();
1084    for (i=baseL_;i<baseL_+numberL_;i++) {
1085      CoinBigIndex j;
1086      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
1087        int iRow = indexRowL[j];
1088        temp[iRow]++;
1089      }
1090    }
1091    for (i=0;i<numberRows_;i++) {
1092      int number = temp[i];
1093      int iPermute = permuteBack[i];
1094      weights[iPermute]=number;
1095    }
1096    delete [] temp;
1097  } else {
1098    int i;
1099    for (i=0;i<numberRows_;i++) {
1100      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
1101      //number = startRowL[i+1]-startRowL[i]+1;
1102      //number = numberInRow[i]+1;
1103      int iPermute = permuteBack[i];
1104      weights[iPermute]=number;
1105    }
1106  }
1107#ifdef CLP_FACTORIZATION_INSTRUMENT
1108  factorization_instrument(8);
1109#endif
1110}
1111#else
1112// This one allows multiple factorizations
1113#if CLP_MULTIPLE_FACTORIZATIONS == 1
1114typedef CoinDenseFactorization CoinOtherFactorization;
1115typedef CoinOslFactorization CoinOtherFactorization;
1116#elif CLP_MULTIPLE_FACTORIZATIONS == 2
1117#include "CoinSimpFactorization.hpp"
1118typedef CoinSimpFactorization CoinOtherFactorization;
1119typedef CoinOslFactorization CoinOtherFactorization;
1120#elif CLP_MULTIPLE_FACTORIZATIONS == 3
1121#include "CoinSimpFactorization.hpp"
1122#define CoinOslFactorization CoinDenseFactorization
1123#elif CLP_MULTIPLE_FACTORIZATIONS == 4
1124#include "CoinSimpFactorization.hpp"
1125#define CoinOslFactorization CoinDenseFactorization
1126//#include "CoinOslFactorization.hpp"
1127#endif
1128
1129//-------------------------------------------------------------------
1130// Default Constructor
1131//-------------------------------------------------------------------
1132ClpFactorization::ClpFactorization () 
1133{
1134#ifndef SLIM_CLP
1135  networkBasis_ = NULL;
1136#endif
1137  //coinFactorizationA_ = NULL;
1138  coinFactorizationA_ = new CoinFactorization() ;
1139  coinFactorizationB_ = NULL;
1140  //coinFactorizationB_ = new CoinOtherFactorization();
1141  forceB_=0;
1142  goOslThreshold_ = -1;
1143  goDenseThreshold_ = -1;
1144  goSmallThreshold_ = -1;
1145}
1146
1147//-------------------------------------------------------------------
1148// Copy constructor
1149//-------------------------------------------------------------------
1150ClpFactorization::ClpFactorization (const ClpFactorization & rhs,
1151                                    int denseIfSmaller) 
1152{
1153#ifdef CLP_FACTORIZATION_INSTRUMENT
1154  factorization_instrument(-1);
1155#endif
1156#ifndef SLIM_CLP
1157  if (rhs.networkBasis_)
1158    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1159  else
1160    networkBasis_=NULL;
1161#endif
1162  forceB_ = rhs.forceB_;
1163  goOslThreshold_ = rhs.goOslThreshold_;
1164  goDenseThreshold_ = rhs.goDenseThreshold_;
1165  goSmallThreshold_ = rhs.goSmallThreshold_;
1166  int goDense = 0;
1167  if (denseIfSmaller>0&&denseIfSmaller<=goDenseThreshold_) {
1168    CoinDenseFactorization * denseR = 
1169      dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
1170    if (!denseR)
1171      goDense=1;
1172  }
1173  if (denseIfSmaller>0&&!rhs.coinFactorizationB_) { 
1174    if (denseIfSmaller<=goDenseThreshold_) 
1175      goDense=1;
1176    else if (denseIfSmaller<=goSmallThreshold_) 
1177      goDense=2;
1178    else if (denseIfSmaller<=goOslThreshold_) 
1179      goDense=3;
1180  } else if (denseIfSmaller<0) {
1181    if (-denseIfSmaller<=goDenseThreshold_) 
1182      goDense=1;
1183    else if (-denseIfSmaller<=goSmallThreshold_) 
1184      goDense=2;
1185    else if (-denseIfSmaller<=goOslThreshold_) 
1186      goDense=3;
1187  }
1188  if (rhs.coinFactorizationA_&&!goDense)
1189    coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
1190  else
1191    coinFactorizationA_=NULL;
1192  if (rhs.coinFactorizationB_&&(denseIfSmaller>=0||!goDense))
1193    coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1194  else
1195    coinFactorizationB_=NULL;
1196  if (goDense) {
1197    delete coinFactorizationB_;
1198    if (goDense==1)
1199      coinFactorizationB_ = new CoinDenseFactorization();
1200    else if (goDense==2)
1201      coinFactorizationB_ = new CoinSimpFactorization();
1202    else
1203      coinFactorizationB_ = new CoinOslFactorization();
1204    if (rhs.coinFactorizationA_) {
1205      coinFactorizationB_->maximumPivots(rhs.coinFactorizationA_->maximumPivots());
1206      coinFactorizationB_->pivotTolerance(rhs.coinFactorizationA_->pivotTolerance());
1207      coinFactorizationB_->zeroTolerance(rhs.coinFactorizationA_->zeroTolerance());
1208    } else {
1209      assert (coinFactorizationB_);
1210      coinFactorizationB_->maximumPivots(rhs.coinFactorizationB_->maximumPivots());
1211      coinFactorizationB_->pivotTolerance(rhs.coinFactorizationB_->pivotTolerance());
1212      coinFactorizationB_->zeroTolerance(rhs.coinFactorizationB_->zeroTolerance());
1213    }
1214  }
1215  assert (!coinFactorizationA_||!coinFactorizationB_);
1216#ifdef CLP_FACTORIZATION_INSTRUMENT
1217  factorization_instrument(1);
1218#endif
1219}
1220
1221ClpFactorization::ClpFactorization (const CoinFactorization & rhs) 
1222{
1223#ifdef CLP_FACTORIZATION_INSTRUMENT
1224  factorization_instrument(-1);
1225#endif
1226#ifndef SLIM_CLP
1227  networkBasis_=NULL;
1228#endif
1229  coinFactorizationA_ = new CoinFactorization(rhs);
1230  coinFactorizationB_=NULL;
1231#ifdef CLP_FACTORIZATION_INSTRUMENT
1232  factorization_instrument(1);
1233#endif
1234  forceB_=0;
1235  goOslThreshold_ = -1;
1236  goDenseThreshold_ = -1;
1237  goSmallThreshold_ = -1;
1238  assert (!coinFactorizationA_||!coinFactorizationB_);
1239}
1240
1241ClpFactorization::ClpFactorization (const CoinOtherFactorization & rhs) 
1242{
1243#ifdef CLP_FACTORIZATION_INSTRUMENT
1244  factorization_instrument(-1);
1245#endif
1246#ifndef SLIM_CLP
1247  networkBasis_=NULL;
1248#endif
1249  coinFactorizationA_ = NULL;
1250  coinFactorizationB_ = rhs.clone();
1251  //coinFactorizationB_ = new CoinOtherFactorization(rhs);
1252  forceB_=0;
1253  goOslThreshold_ = -1;
1254  goDenseThreshold_ = -1;
1255  goSmallThreshold_ = -1;
1256#ifdef CLP_FACTORIZATION_INSTRUMENT
1257  factorization_instrument(1);
1258#endif
1259  assert (!coinFactorizationA_||!coinFactorizationB_);
1260}
1261
1262//-------------------------------------------------------------------
1263// Destructor
1264//-------------------------------------------------------------------
1265ClpFactorization::~ClpFactorization () 
1266{
1267#ifndef SLIM_CLP
1268  delete networkBasis_;
1269#endif
1270  delete coinFactorizationA_;
1271  delete coinFactorizationB_;
1272}
1273
1274//----------------------------------------------------------------
1275// Assignment operator
1276//-------------------------------------------------------------------
1277ClpFactorization &
1278ClpFactorization::operator=(const ClpFactorization& rhs)
1279{
1280#ifdef CLP_FACTORIZATION_INSTRUMENT
1281  factorization_instrument(-1);
1282#endif
1283  if (this != &rhs) {
1284#ifndef SLIM_CLP
1285    delete networkBasis_;
1286    if (rhs.networkBasis_)
1287      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1288    else
1289      networkBasis_=NULL;
1290#endif
1291    forceB_ = rhs.forceB_;
1292    goOslThreshold_ = rhs.goOslThreshold_;
1293    goDenseThreshold_ = rhs.goDenseThreshold_;
1294    goSmallThreshold_ = rhs.goSmallThreshold_;
1295    if (rhs.coinFactorizationA_) {
1296      if (coinFactorizationA_)
1297        *coinFactorizationA_ = *(rhs.coinFactorizationA_);
1298      else
1299        coinFactorizationA_ = new CoinFactorization(*rhs.coinFactorizationA_);
1300    } else {
1301      delete coinFactorizationA_;
1302      coinFactorizationA_=NULL;
1303    }
1304   
1305    if (rhs.coinFactorizationB_) {
1306      if (coinFactorizationB_) {
1307        CoinDenseFactorization * denseR = dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
1308        CoinDenseFactorization * dense = dynamic_cast<CoinDenseFactorization *>(coinFactorizationB_);
1309        CoinOslFactorization * oslR = dynamic_cast<CoinOslFactorization *>(rhs.coinFactorizationB_);
1310        CoinOslFactorization * osl = dynamic_cast<CoinOslFactorization *>(coinFactorizationB_);
1311        CoinSimpFactorization * simpR = dynamic_cast<CoinSimpFactorization *>(rhs.coinFactorizationB_);
1312        CoinSimpFactorization * simp = dynamic_cast<CoinSimpFactorization *>(coinFactorizationB_);
1313        if (dense&&denseR) {
1314          *dense=*denseR;
1315        } else if (osl&&oslR) {
1316          *osl=*oslR;
1317        } else if (simp&&simpR) {
1318          *simp=*simpR;
1319        } else {
1320          delete coinFactorizationB_;
1321          coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1322        }
1323      } else {
1324        coinFactorizationB_ = rhs.coinFactorizationB_->clone();
1325      }
1326    } else {
1327      delete coinFactorizationB_;
1328      coinFactorizationB_=NULL;
1329    }
1330  }
1331#ifdef CLP_FACTORIZATION_INSTRUMENT
1332  factorization_instrument(1);
1333#endif
1334  assert (!coinFactorizationA_||!coinFactorizationB_);
1335  return *this;
1336}
1337// Go over to dense code
1338void 
1339ClpFactorization::goDenseOrSmall(int numberRows) 
1340{
1341  if (!forceB_) {
1342    if (numberRows<=goDenseThreshold_) {
1343      delete coinFactorizationA_;
1344      delete coinFactorizationB_;
1345      coinFactorizationA_ = NULL;
1346      coinFactorizationB_ = new CoinDenseFactorization();
1347      //printf("going dense\n");
1348    } else if (numberRows<=goSmallThreshold_) {
1349      delete coinFactorizationA_;
1350      delete coinFactorizationB_;
1351      coinFactorizationA_ = NULL;
1352      coinFactorizationB_ = new CoinSimpFactorization();
1353      //printf("going small\n");
1354    } else if (numberRows<=goOslThreshold_) {
1355      delete coinFactorizationA_;
1356      delete coinFactorizationB_;
1357      coinFactorizationA_ = NULL;
1358      coinFactorizationB_ = new CoinOslFactorization();
1359      //printf("going small\n");
1360    }
1361  }
1362  assert (!coinFactorizationA_||!coinFactorizationB_);
1363}
1364// If nonzero force use of 1,dense 2,small 3,osl
1365void 
1366ClpFactorization::forceOtherFactorization(int which)
1367{
1368  delete coinFactorizationB_;
1369  forceB_=0;
1370  coinFactorizationB_ = NULL;
1371  if (which>0&&which<4) {
1372    delete coinFactorizationA_;
1373    coinFactorizationA_ = NULL;
1374    forceB_=which;
1375    switch (which) {
1376    case 1:
1377      coinFactorizationB_ = new CoinDenseFactorization();
1378      goDenseThreshold_=COIN_INT_MAX;
1379    break;
1380    case 2:
1381      coinFactorizationB_ = new CoinSimpFactorization();
1382      goSmallThreshold_=COIN_INT_MAX;
1383    break;
1384    case 3:
1385      coinFactorizationB_ = new CoinOslFactorization();
1386      goOslThreshold_=COIN_INT_MAX;
1387    break;
1388    }
1389  } else if (!coinFactorizationA_) {
1390    coinFactorizationA_ = new CoinFactorization();
1391  }
1392}
1393int 
1394ClpFactorization::factorize ( ClpSimplex * model,
1395                              int solveType, bool valuesPass)
1396{
1397  //if ((model->specialOptions()&16384))
1398  //printf("factor at %d iterations\n",model->numberIterations());
1399  ClpMatrixBase * matrix = model->clpMatrix(); 
1400  int numberRows = model->numberRows();
1401  int numberColumns = model->numberColumns();
1402  if (!numberRows)
1403    return 0;
1404#ifdef CLP_FACTORIZATION_INSTRUMENT
1405  factorization_instrument(-1);
1406#endif
1407  bool anyChanged=false;
1408  if (coinFactorizationB_) {
1409    coinFactorizationB_->setStatus(-99);
1410    int * pivotVariable = model->pivotVariable();
1411    //returns 0 -okay, -1 singular, -2 too many in basis */
1412    // allow dense
1413    int solveMode=2;
1414    if (model->numberIterations())
1415      solveMode+=8;
1416    if (valuesPass)
1417      solveMode+= 4;
1418    coinFactorizationB_->setSolveMode(solveMode);
1419    while (status()<-98) {
1420     
1421      int i;
1422      int numberBasic=0;
1423      int numberRowBasic;
1424      // Move pivot variables across if they look good
1425      int * pivotTemp = model->rowArray(0)->getIndices();
1426      assert (!model->rowArray(0)->getNumElements());
1427      if (!matrix->rhsOffset(model)) {
1428        // Seems to prefer things in order so quickest
1429        // way is to go though like this
1430        for (i=0;i<numberRows;i++) {
1431          if (model->getRowStatus(i) == ClpSimplex::basic) 
1432            pivotTemp[numberBasic++]=i;
1433        }
1434        numberRowBasic=numberBasic;
1435        /* Put column basic variables into pivotVariable
1436           This is done by ClpMatrixBase to allow override for gub
1437        */
1438        matrix->generalExpanded(model,0,numberBasic);
1439      } else {
1440        // Long matrix - do a different way
1441        bool fullSearch=false;
1442        for (i=0;i<numberRows;i++) {
1443          int iPivot = pivotVariable[i];
1444          if (iPivot>=numberColumns) {
1445            pivotTemp[numberBasic++]=iPivot-numberColumns;
1446          }
1447        }
1448        numberRowBasic=numberBasic;
1449        for (i=0;i<numberRows;i++) {
1450          int iPivot = pivotVariable[i];
1451          if (iPivot<numberColumns) {
1452            if (iPivot>=0) {
1453              pivotTemp[numberBasic++]=iPivot;
1454            } else {
1455              // not full basis
1456              fullSearch=true;
1457              break;
1458            }
1459          }
1460        }
1461        if (fullSearch) {
1462          // do slow way
1463          numberBasic=0;
1464          for (i=0;i<numberRows;i++) {
1465            if (model->getRowStatus(i) == ClpSimplex::basic) 
1466              pivotTemp[numberBasic++]=i;
1467          }
1468          numberRowBasic=numberBasic;
1469          /* Put column basic variables into pivotVariable
1470             This is done by ClpMatrixBase to allow override for gub
1471          */
1472          matrix->generalExpanded(model,0,numberBasic);
1473        }
1474      }
1475      if (numberBasic>model->maximumBasic()) {
1476        // Take out some
1477        numberBasic=numberRowBasic;
1478        for (int i=0;i<numberColumns;i++) {
1479          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1480            if (numberBasic<numberRows)
1481              numberBasic++;
1482            else
1483              model->setColumnStatus(i,ClpSimplex::superBasic);
1484          }
1485        }
1486        numberBasic=numberRowBasic;
1487        matrix->generalExpanded(model,0,numberBasic);
1488      } else if (numberBasic<numberRows) {
1489        // add in some
1490        int needed = numberRows-numberBasic;
1491        // move up columns
1492        for (i=numberBasic-1;i>=numberRowBasic;i--)
1493          pivotTemp[i+needed]=pivotTemp[i];
1494        numberRowBasic=0;
1495        numberBasic=numberRows;
1496        for (i=0;i<numberRows;i++) {
1497          if (model->getRowStatus(i) == ClpSimplex::basic) {
1498            pivotTemp[numberRowBasic++]=i;
1499          } else if (needed) {
1500            needed--;
1501            model->setRowStatus(i,ClpSimplex::basic);
1502            pivotTemp[numberRowBasic++]=i;
1503          }
1504        }
1505      }
1506      CoinBigIndex numberElements=numberRowBasic;
1507     
1508      // compute how much in basis
1509      // can change for gub
1510      int numberColumnBasic = numberBasic-numberRowBasic;
1511     
1512      numberElements +=matrix->countBasis(model,
1513                                          pivotTemp+numberRowBasic, 
1514                                          numberRowBasic,
1515                                          numberColumnBasic);
1516      // Not needed for dense
1517      numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1518      int numberIterations=model->numberIterations();
1519      coinFactorizationB_->setUsefulInformation(&numberIterations,0);
1520      coinFactorizationB_->getAreas ( numberRows,
1521                                      numberRowBasic+numberColumnBasic, numberElements,
1522                                      2 * numberElements );
1523      // Fill in counts so we can skip part of preProcess
1524      // This is NOT needed for dense but would be needed for later versions
1525      CoinFactorizationDouble * elementU;
1526      int * indexRowU;
1527      CoinBigIndex * startColumnU;
1528      int * numberInRow;
1529      int * numberInColumn;
1530      elementU = coinFactorizationB_->elements();
1531      indexRowU = coinFactorizationB_->indices();
1532      startColumnU = coinFactorizationB_->starts();
1533#ifndef COIN_FAST_CODE
1534      double slackValue;
1535      slackValue = coinFactorizationB_->slackValue();
1536#else
1537#define slackValue -1.0
1538#endif
1539      numberInRow = coinFactorizationB_->numberInRow();
1540      numberInColumn = coinFactorizationB_->numberInColumn();
1541      CoinZeroN ( numberInRow, numberRows  );
1542      CoinZeroN ( numberInColumn, numberRows );
1543      for (i=0;i<numberRowBasic;i++) {
1544        int iRow = pivotTemp[i];
1545        // Change pivotTemp to correct sequence
1546        pivotTemp[i]=iRow+numberColumns;
1547        indexRowU[i]=iRow;
1548        startColumnU[i]=i;
1549        elementU[i]=slackValue;
1550        numberInRow[iRow]=1;
1551        numberInColumn[i]=1;
1552      }
1553      startColumnU[numberRowBasic]=numberRowBasic;
1554      // can change for gub so redo
1555      numberColumnBasic = numberBasic-numberRowBasic;
1556      matrix->fillBasis(model, 
1557                        pivotTemp+numberRowBasic, 
1558                        numberColumnBasic,
1559                        indexRowU, 
1560                        startColumnU+numberRowBasic,
1561                        numberInRow,
1562                        numberInColumn+numberRowBasic,
1563                        elementU);
1564      // recompute number basic
1565      numberBasic = numberRowBasic+numberColumnBasic;
1566      for (i=numberBasic;i<numberRows;i++)
1567        pivotTemp[i]=-1; // mark not there
1568      if (numberBasic) 
1569        numberElements = startColumnU[numberBasic-1]
1570          +numberInColumn[numberBasic-1];
1571      else
1572        numberElements=0;
1573      coinFactorizationB_->preProcess ( );
1574      coinFactorizationB_->factor (  );
1575      if (coinFactorizationB_->status() == -1 &&
1576          (coinFactorizationB_->solveMode()%3)!=0) {
1577        int solveMode=coinFactorizationB_->solveMode();
1578        solveMode -= solveMode%3; // so bottom will be 0
1579        coinFactorizationB_->setSolveMode(solveMode);
1580        coinFactorizationB_->setStatus(-99);
1581      }
1582      if (coinFactorizationB_->status() == -99)
1583        continue;
1584      // If we get here status is 0 or -1
1585      if (coinFactorizationB_->status() == 0&&numberBasic==numberRows) {
1586        coinFactorizationB_->postProcess(pivotTemp,pivotVariable);
1587      } else if (solveType==0||solveType==2/*||solveType==1*/) {
1588        // Change pivotTemp to be correct list
1589        anyChanged=true;
1590        coinFactorizationB_->makeNonSingular(pivotTemp,numberColumns);
1591        double * columnLower = model->lowerRegion();
1592        double * columnUpper = model->upperRegion();
1593        double * columnActivity = model->solutionRegion();
1594        double * rowLower = model->lowerRegion(0);
1595        double * rowUpper = model->upperRegion(0);
1596        double * rowActivity = model->solutionRegion(0);
1597        //redo basis - first take ALL out
1598        int iColumn;
1599        double largeValue = model->largeValue();
1600        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1601          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
1602            // take out
1603            if (!valuesPass) {
1604              double lower=columnLower[iColumn];
1605              double upper=columnUpper[iColumn];
1606              double value=columnActivity[iColumn];
1607              if (lower>-largeValue||upper<largeValue) {
1608                if (fabs(value-lower)<fabs(value-upper)) {
1609                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
1610                  columnActivity[iColumn]=lower;
1611                } else {
1612                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
1613                  columnActivity[iColumn]=upper;
1614                }
1615              } else {
1616                model->setColumnStatus(iColumn,ClpSimplex::isFree);
1617              }
1618            } else {
1619              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
1620            }
1621          }
1622        }
1623        int iRow;
1624        for (iRow=0;iRow<numberRows;iRow++) {
1625          if (model->getRowStatus(iRow)==ClpSimplex::basic) {
1626            // take out
1627            if (!valuesPass) {
1628              double lower=columnLower[iRow];
1629              double upper=columnUpper[iRow];
1630              double value=columnActivity[iRow];
1631              if (lower>-largeValue||upper<largeValue) {
1632                if (fabs(value-lower)<fabs(value-upper)) {
1633                  model->setRowStatus(iRow,ClpSimplex::atLowerBound);
1634                  columnActivity[iRow]=lower;
1635                } else {
1636                  model->setRowStatus(iRow,ClpSimplex::atUpperBound);
1637                  columnActivity[iRow]=upper;
1638                }
1639              } else {
1640                model->setRowStatus(iRow,ClpSimplex::isFree);
1641              }
1642            } else {
1643              model->setRowStatus(iRow,ClpSimplex::superBasic);
1644            }
1645          }
1646        }
1647        for (iRow=0;iRow<numberRows;iRow++) {
1648          int iSequence=pivotTemp[iRow];
1649          assert (iSequence>=0);
1650          // basic
1651          model->setColumnStatus(iSequence,ClpSimplex::basic);
1652        }
1653        // signal repeat
1654        coinFactorizationB_->setStatus(-99);
1655        // set fixed if they are
1656        for (iRow=0;iRow<numberRows;iRow++) {
1657          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
1658            if (rowLower[iRow]==rowUpper[iRow]) {
1659              rowActivity[iRow]=rowLower[iRow];
1660              model->setRowStatus(iRow,ClpSimplex::isFixed);
1661            }
1662          }
1663        }
1664        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1665          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
1666            if (columnLower[iColumn]==columnUpper[iColumn]) {
1667              columnActivity[iColumn]=columnLower[iColumn];
1668              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
1669            }
1670          }
1671        }
1672      } 
1673    }
1674#ifdef CLP_DEBUG
1675    // check basic
1676    CoinIndexedVector region1(2*numberRows);
1677    CoinIndexedVector region2B(2*numberRows);
1678    int iPivot;
1679    double * arrayB = region2B.denseVector();
1680    int i;
1681    for (iPivot=0;iPivot<numberRows;iPivot++) {
1682      int iSequence = pivotVariable[iPivot];
1683      model->unpack(&region2B,iSequence);
1684      coinFactorizationB_->updateColumn(&region1,&region2B);
1685      if (fabs(arrayB[iPivot]-1.0)<1.0e-4) {
1686        // OK?
1687        arrayB[iPivot]=0.0;
1688      } else {
1689        assert (fabs(arrayB[iPivot])<1.0e-4);
1690        for (i=0;i<numberRows;i++) {
1691          if (fabs(arrayB[i]-1.0)<1.0e-4)
1692            break;
1693        }
1694        assert (i<numberRows);
1695        printf("variable on row %d landed up on row %d\n",iPivot,i);
1696        arrayB[i]=0.0;
1697      }
1698      for (i=0;i<numberRows;i++)
1699        assert (fabs(arrayB[i])<1.0e-4);
1700      region2B.clear();
1701    }
1702#endif
1703#ifdef CLP_FACTORIZATION_INSTRUMENT
1704    factorization_instrument(2);
1705#endif
1706    if ( anyChanged&&model->algorithm()<0&&solveType>0) {
1707      double dummyCost;
1708      static_cast<ClpSimplexDual *> (model)->changeBounds(3,
1709                                                          NULL,dummyCost);
1710    }
1711    return coinFactorizationB_->status();
1712  }
1713  // If too many compressions increase area
1714  if (coinFactorizationA_->pivots()>1&&coinFactorizationA_->numberCompressions()*10 > coinFactorizationA_->pivots()+10) {
1715    coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor()* 1.1);
1716  }
1717  //int numberPivots=coinFactorizationA_->pivots();
1718#if 0
1719  if (model->algorithm()>0)
1720    numberSave=-1;
1721#endif
1722#ifndef SLIM_CLP
1723  if (!networkBasis_||doCheck) {
1724#endif
1725    coinFactorizationA_->setStatus(-99);
1726    int * pivotVariable = model->pivotVariable();
1727    int nTimesRound=0;
1728    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
1729    while (coinFactorizationA_->status()<-98) {
1730      nTimesRound++;
1731     
1732      int i;
1733      int numberBasic=0;
1734      int numberRowBasic;
1735      // Move pivot variables across if they look good
1736      int * pivotTemp = model->rowArray(0)->getIndices();
1737      assert (!model->rowArray(0)->getNumElements());
1738      if (!matrix->rhsOffset(model)) {
1739#if 0
1740        if (numberSave>0) {
1741          int nStill=0;
1742          int nAtBound=0;
1743          int nZeroDual=0;
1744          CoinIndexedVector * array = model->rowArray(3);
1745          CoinIndexedVector * objArray = model->columnArray(1);
1746          array->clear();
1747          objArray->clear();
1748          double * cost = model->costRegion();
1749          double tolerance = model->primalTolerance();
1750          double offset=0.0;
1751          for (i=0;i<numberRows;i++) {
1752            int iPivot = pivotVariable[i];
1753            if (iPivot<numberColumns&&isDense(iPivot)) {
1754              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
1755                nStill++;
1756                double value=model->solutionRegion()[iPivot];
1757                double dual = model->dualRowSolution()[i];
1758                double lower=model->lowerRegion()[iPivot];
1759                double upper=model->upperRegion()[iPivot];
1760                ClpSimplex::Status status;
1761                if (fabs(value-lower)<tolerance) {
1762                  status=ClpSimplex::atLowerBound;
1763                  nAtBound++;
1764                } else if (fabs(value-upper)<tolerance) {
1765                  nAtBound++;
1766                  status=ClpSimplex::atUpperBound;
1767                } else if (value>lower&&value<upper) {
1768                  status=ClpSimplex::superBasic;
1769                } else {
1770                  status=ClpSimplex::basic;
1771                }
1772                if (status!=ClpSimplex::basic) {
1773                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
1774                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
1775                    model->setRowStatus(i,ClpSimplex::basic);
1776                    pivotVariable[i]=i+numberColumns;
1777                    model->dualRowSolution()[i]=0.0;
1778                    model->djRegion(0)[i]=0.0;
1779                    array->add(i,dual);
1780                    offset += dual*model->solutionRegion(0)[i];
1781                  }
1782                }
1783                if (fabs(dual)<1.0e-5)
1784                  nZeroDual++;
1785              }
1786            }
1787          }
1788          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
1789                 numberSave,nStill,nAtBound,nZeroDual,offset);
1790          if (array->getNumElements()) {
1791            // modify costs
1792            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
1793                                               objArray);
1794            array->clear();
1795            int n=objArray->getNumElements();
1796            int * indices = objArray->getIndices();
1797            double * elements = objArray->denseVector();
1798            for (i=0;i<n;i++) {
1799              int iColumn = indices[i];
1800              cost[iColumn] -= elements[iColumn];
1801              elements[iColumn]=0.0;
1802            }
1803            objArray->setNumElements(0);
1804          }
1805        }
1806#endif
1807        // Seems to prefer things in order so quickest
1808        // way is to go though like this
1809        for (i=0;i<numberRows;i++) {
1810          if (model->getRowStatus(i) == ClpSimplex::basic) 
1811            pivotTemp[numberBasic++]=i;
1812        }
1813        numberRowBasic=numberBasic;
1814        /* Put column basic variables into pivotVariable
1815           This is done by ClpMatrixBase to allow override for gub
1816        */
1817        matrix->generalExpanded(model,0,numberBasic);
1818      } else {
1819        // Long matrix - do a different way
1820        bool fullSearch=false;
1821        for (i=0;i<numberRows;i++) {
1822          int iPivot = pivotVariable[i];
1823          if (iPivot>=numberColumns) {
1824            pivotTemp[numberBasic++]=iPivot-numberColumns;
1825          }
1826        }
1827        numberRowBasic=numberBasic;
1828        for (i=0;i<numberRows;i++) {
1829          int iPivot = pivotVariable[i];
1830          if (iPivot<numberColumns) {
1831            if (iPivot>=0) {
1832              pivotTemp[numberBasic++]=iPivot;
1833            } else {
1834              // not full basis
1835              fullSearch=true;
1836              break;
1837            }
1838          }
1839        }
1840        if (fullSearch) {
1841          // do slow way
1842          numberBasic=0;
1843          for (i=0;i<numberRows;i++) {
1844            if (model->getRowStatus(i) == ClpSimplex::basic) 
1845              pivotTemp[numberBasic++]=i;
1846          }
1847          numberRowBasic=numberBasic;
1848          /* Put column basic variables into pivotVariable
1849             This is done by ClpMatrixBase to allow override for gub
1850          */
1851          matrix->generalExpanded(model,0,numberBasic);
1852        }
1853      }
1854      if (numberBasic>model->maximumBasic()) {
1855#if 0 // ndef NDEBUG
1856        printf("%d basic - should only be %d\n",
1857               numberBasic,numberRows);
1858#endif
1859        // Take out some
1860        numberBasic=numberRowBasic;
1861        for (int i=0;i<numberColumns;i++) {
1862          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1863            if (numberBasic<numberRows)
1864              numberBasic++;
1865            else
1866              model->setColumnStatus(i,ClpSimplex::superBasic);
1867          }
1868        }
1869        numberBasic=numberRowBasic;
1870        matrix->generalExpanded(model,0,numberBasic);
1871      }
1872#ifndef SLIM_CLP
1873      // see if matrix a network
1874#ifndef NO_RTTI
1875      ClpNetworkMatrix* networkMatrix =
1876        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
1877#else
1878      ClpNetworkMatrix* networkMatrix = NULL;
1879      if (model->clpMatrix()->type()==11)
1880        networkMatrix = 
1881        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
1882#endif
1883      // If network - still allow ordinary factorization first time for laziness
1884      if (networkMatrix)
1885        coinFactorizationA_->setBiasLU(0); // All to U if network
1886      //int saveMaximumPivots = maximumPivots();
1887      delete networkBasis_;
1888      networkBasis_ = NULL;
1889      if (networkMatrix&&!doCheck)
1890        maximumPivots(1);
1891#endif
1892      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
1893      while (coinFactorizationA_->status()==-99) {
1894        // maybe for speed will be better to leave as many regions as possible
1895        coinFactorizationA_->gutsOfDestructor();
1896        coinFactorizationA_->gutsOfInitialize(2);
1897        CoinBigIndex numberElements=numberRowBasic;
1898
1899        // compute how much in basis
1900
1901        int i;
1902        // can change for gub
1903        int numberColumnBasic = numberBasic-numberRowBasic;
1904
1905        numberElements +=matrix->countBasis(model,
1906                                           pivotTemp+numberRowBasic, 
1907                                           numberRowBasic,
1908                                            numberColumnBasic);
1909        // and recompute as network side say different
1910        if (model->numberIterations())
1911          numberRowBasic = numberBasic - numberColumnBasic;
1912        numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1913        coinFactorizationA_->getAreas ( numberRows,
1914                   numberRowBasic+numberColumnBasic, numberElements,
1915                   2 * numberElements );
1916        //fill
1917        // Fill in counts so we can skip part of preProcess
1918        int * numberInRow = coinFactorizationA_->numberInRow();
1919        int * numberInColumn = coinFactorizationA_->numberInColumn();
1920        CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
1921        CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
1922        CoinFactorizationDouble * elementU = coinFactorizationA_->elementU();
1923        int * indexRowU = coinFactorizationA_->indexRowU();
1924        CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
1925#ifndef COIN_FAST_CODE
1926        double slackValue = coinFactorizationA_->slackValue();
1927#endif
1928        for (i=0;i<numberRowBasic;i++) {
1929          int iRow = pivotTemp[i];
1930          indexRowU[i]=iRow;
1931          startColumnU[i]=i;
1932          elementU[i]=slackValue;
1933          numberInRow[iRow]=1;
1934          numberInColumn[i]=1;
1935        }
1936        startColumnU[numberRowBasic]=numberRowBasic;
1937        // can change for gub so redo
1938        numberColumnBasic = numberBasic-numberRowBasic;
1939        matrix->fillBasis(model, 
1940                          pivotTemp+numberRowBasic, 
1941                          numberColumnBasic,
1942                          indexRowU, 
1943                          startColumnU+numberRowBasic,
1944                          numberInRow,
1945                          numberInColumn+numberRowBasic,
1946                          elementU);
1947#if 0
1948        {
1949          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
1950          for (int i=0;i<numberElements;i++)
1951            printf("row %d col %d value %g\n",indexRowU[i],indexColumnU_[i],
1952                   elementU[i]);
1953        }
1954#endif
1955        // recompute number basic
1956        numberBasic = numberRowBasic+numberColumnBasic;
1957        if (numberBasic) 
1958          numberElements = startColumnU[numberBasic-1]
1959            +numberInColumn[numberBasic-1];
1960        else
1961          numberElements=0;
1962        coinFactorizationA_->setNumberElementsU(numberElements);
1963        //saveFactorization("dump.d");
1964        if (coinFactorizationA_->biasLU()>=3||coinFactorizationA_->numberRows()!=coinFactorizationA_->numberColumns())
1965          coinFactorizationA_->preProcess ( 2 );
1966        else
1967          coinFactorizationA_->preProcess ( 3 ); // no row copy
1968        coinFactorizationA_->factor (  );
1969        if (coinFactorizationA_->status()==-99) {
1970          // get more memory
1971          coinFactorizationA_->areaFactor(2.0*coinFactorizationA_->areaFactor());
1972        } else if (coinFactorizationA_->status()==-1&&
1973                   (model->numberIterations()==0||nTimesRound>2)&&
1974                   coinFactorizationA_->denseThreshold()) {
1975          // Round again without dense
1976          coinFactorizationA_->setDenseThreshold(0);
1977          coinFactorizationA_->setStatus(-99);
1978        }
1979      }
1980      // If we get here status is 0 or -1
1981      if (coinFactorizationA_->status() == 0) {
1982        // We may need to tamper with order and redo - e.g. network with side
1983        int useNumberRows = numberRows;
1984        // **** we will also need to add test in dual steepest to do
1985        // as we do for network
1986        matrix->generalExpanded(model,12,useNumberRows);
1987        const int * permuteBack = coinFactorizationA_->permuteBack();
1988        const int * back = coinFactorizationA_->pivotColumnBack();
1989        //int * pivotTemp = pivotColumn_.array();
1990        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
1991#ifndef NDEBUG
1992        CoinFillN(pivotVariable,numberRows,-1);
1993#endif
1994        // Redo pivot order
1995        for (i=0;i<numberRowBasic;i++) {
1996          int k = pivotTemp[i];
1997          // so rowIsBasic[k] would be permuteBack[back[i]]
1998          int j=permuteBack[back[i]];
1999          assert (pivotVariable[j]==-1);
2000          pivotVariable[j]=k+numberColumns;
2001        }
2002        for (;i<useNumberRows;i++) {
2003          int k = pivotTemp[i];
2004          // so rowIsBasic[k] would be permuteBack[back[i]]
2005          int j=permuteBack[back[i]];
2006          assert (pivotVariable[j]==-1);
2007          pivotVariable[j]=k;
2008        }
2009#if 0
2010        if (numberSave>=0) {
2011          numberSave=numberDense_;
2012          memset(saveList,0,((coinFactorizationA_->numberRows()+31)>>5)*sizeof(int));
2013          for (i=coinFactorizationA_->numberRows()-numberSave;i<coinFactorizationA_->numberRows();i++) {
2014            int k=pivotTemp[pivotColumn_.array()[i]];
2015            setDense(k);
2016          }
2017        }
2018#endif
2019        // Set up permutation vector
2020        // these arrays start off as copies of permute
2021        // (and we could use permute_ instead of pivotColumn (not back though))
2022        ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn()  );
2023        ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack()  );
2024#ifndef SLIM_CLP
2025        if (networkMatrix) {
2026          maximumPivots(CoinMax(2000,maximumPivots()));
2027          // redo arrays
2028          for (int iRow=0;iRow<4;iRow++) {
2029            int length =model->numberRows()+maximumPivots();
2030            if (iRow==3||model->objectiveAsObject()->type()>1)
2031              length += model->numberColumns();
2032            model->rowArray(iRow)->reserve(length);
2033          }
2034          // create network factorization
2035          if (doCheck)
2036            delete networkBasis_; // temp
2037          networkBasis_ = new ClpNetworkBasis(model,coinFactorizationA_->numberRows(),
2038                                              coinFactorizationA_->pivotRegion(),
2039                                              coinFactorizationA_->permuteBack(),
2040                                              coinFactorizationA_->startColumnU(),
2041                                              coinFactorizationA_->numberInColumn(),
2042                                              coinFactorizationA_->indexRowU(),
2043                                              coinFactorizationA_->elementU());
2044          // kill off arrays in ordinary factorization
2045          if (!doCheck) {
2046            coinFactorizationA_->gutsOfDestructor();
2047            // but make sure coinFactorizationA_->numberRows() set
2048            coinFactorizationA_->setNumberRows(model->numberRows());
2049            coinFactorizationA_->setStatus(0);
2050#if 0
2051            // but put back permute arrays so odd things will work
2052            int numberRows = model->numberRows();
2053            pivotColumnBack_ = new int [numberRows];
2054            permute_ = new int [numberRows];
2055            int i;
2056            for (i=0;i<numberRows;i++) {
2057              pivotColumnBack_[i]=i;
2058              permute_[i]=i;
2059            }
2060#endif
2061          }
2062        } else {
2063#endif
2064          // See if worth going sparse and when
2065          coinFactorizationA_->checkSparse();
2066#ifndef SLIM_CLP
2067        }
2068#endif
2069      } else if (coinFactorizationA_->status() == -1&&(solveType==0||solveType==2)) {
2070        // This needs redoing as it was merged coding - does not need array
2071        int numberTotal = numberRows+numberColumns;
2072        int * isBasic = new int [numberTotal];
2073        int * rowIsBasic = isBasic+numberColumns;
2074        int * columnIsBasic = isBasic;
2075        for (i=0;i<numberTotal;i++) 
2076          isBasic[i]=-1;
2077        for (i=0;i<numberRowBasic;i++) {
2078          int iRow = pivotTemp[i];
2079          rowIsBasic[iRow]=1;
2080        }
2081        for (;i<numberBasic;i++) {
2082          int iColumn = pivotTemp[i];
2083          columnIsBasic[iColumn]=1;
2084        }
2085        numberBasic=0;
2086        for (i=0;i<numberRows;i++) 
2087          pivotVariable[i]=-1;
2088        // mark as basic or non basic
2089        const int * pivotColumn = coinFactorizationA_->pivotColumn();
2090        for (i=0;i<numberRows;i++) {
2091          if (rowIsBasic[i]>=0) {
2092            if (pivotColumn[numberBasic]>=0) {
2093              rowIsBasic[i]=pivotColumn[numberBasic];
2094            } else {
2095              rowIsBasic[i]=-1;
2096              model->setRowStatus(i,ClpSimplex::superBasic);
2097            }
2098            numberBasic++;
2099          }
2100        }
2101        for (i=0;i<numberColumns;i++) {
2102          if (columnIsBasic[i]>=0) {
2103            if (pivotColumn[numberBasic]>=0) 
2104              columnIsBasic[i]=pivotColumn[numberBasic];
2105            else
2106              columnIsBasic[i]=-1;
2107            numberBasic++;
2108          }
2109        }
2110        // leave pivotVariable in useful form for cleaning basis
2111        int * pivotVariable = model->pivotVariable();
2112        for (i=0;i<numberRows;i++) {
2113          pivotVariable[i]=-1;
2114        }
2115
2116        for (i=0;i<numberRows;i++) {
2117          if (model->getRowStatus(i) == ClpSimplex::basic) {
2118            int iPivot = rowIsBasic[i];
2119            if (iPivot>=0) 
2120              pivotVariable[iPivot]=i+numberColumns;
2121          }
2122        }
2123        for (i=0;i<numberColumns;i++) {
2124          if (model->getColumnStatus(i) == ClpSimplex::basic) {
2125            int iPivot = columnIsBasic[i];
2126            if (iPivot>=0) 
2127              pivotVariable[iPivot]=i;
2128          }
2129        }
2130        delete [] isBasic;
2131        double * columnLower = model->lowerRegion();
2132        double * columnUpper = model->upperRegion();
2133        double * columnActivity = model->solutionRegion();
2134        double * rowLower = model->lowerRegion(0);
2135        double * rowUpper = model->upperRegion(0);
2136        double * rowActivity = model->solutionRegion(0);
2137        //redo basis - first take ALL columns out
2138        int iColumn;
2139        double largeValue = model->largeValue();
2140        for (iColumn=0;iColumn<numberColumns;iColumn++) {
2141          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
2142            // take out
2143            if (!valuesPass) {
2144              double lower=columnLower[iColumn];
2145              double upper=columnUpper[iColumn];
2146              double value=columnActivity[iColumn];
2147              if (lower>-largeValue||upper<largeValue) {
2148                if (fabs(value-lower)<fabs(value-upper)) {
2149                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
2150                  columnActivity[iColumn]=lower;
2151                } else {
2152                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
2153                  columnActivity[iColumn]=upper;
2154                }
2155              } else {
2156                model->setColumnStatus(iColumn,ClpSimplex::isFree);
2157              }
2158            } else {
2159              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
2160            }
2161          }
2162        }
2163        int iRow;
2164        for (iRow=0;iRow<numberRows;iRow++) {
2165          int iSequence=pivotVariable[iRow];
2166          if (iSequence>=0) {
2167            // basic
2168            if (iSequence>=numberColumns) {
2169              // slack in - leave
2170              //assert (iSequence-numberColumns==iRow);
2171            } else {
2172              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
2173              // put back structural
2174              model->setColumnStatus(iSequence,ClpSimplex::basic);
2175            }
2176          } else {
2177            // put in slack
2178            model->setRowStatus(iRow,ClpSimplex::basic);
2179          }
2180        }
2181        // Put back any key variables for gub
2182        int dummy;
2183        matrix->generalExpanded(model,1,dummy);
2184        // signal repeat
2185        coinFactorizationA_->setStatus(-99);
2186        // set fixed if they are
2187        for (iRow=0;iRow<numberRows;iRow++) {
2188          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
2189            if (rowLower[iRow]==rowUpper[iRow]) {
2190              rowActivity[iRow]=rowLower[iRow];
2191              model->setRowStatus(iRow,ClpSimplex::isFixed);
2192            }
2193          }
2194        }
2195        for (iColumn=0;iColumn<numberColumns;iColumn++) {
2196          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
2197            if (columnLower[iColumn]==columnUpper[iColumn]) {
2198              columnActivity[iColumn]=columnLower[iColumn];
2199              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
2200            }
2201          }
2202        }
2203      } 
2204    }
2205#ifndef SLIM_CLP
2206  } else {
2207    // network - fake factorization - do nothing
2208    coinFactorizationA_->setStatus(0);
2209    coinFactorizationA_->setPivots(0);
2210  }
2211#endif
2212#ifndef SLIM_CLP
2213  if (!coinFactorizationA_->status()) {
2214    // take out part if quadratic
2215    if (model->algorithm()==2) {
2216      ClpObjective * obj = model->objectiveAsObject();
2217#ifndef NDEBUG
2218      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2219      assert (quadraticObj);
2220#else
2221      ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
2222#endif
2223      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2224      int numberXColumns = quadratic->getNumCols();
2225      assert (numberXColumns<numberColumns);
2226      int base = numberColumns-numberXColumns;
2227      int * which = new int [numberXColumns];
2228      int * pivotVariable = model->pivotVariable();
2229      int * permute = pivotColumn();
2230      int i;
2231      int n=0;
2232      for (i=0;i<numberRows;i++) {
2233        int iSj = pivotVariable[i]-base;
2234        if (iSj>=0&&iSj<numberXColumns) 
2235          which[n++]=permute[i];
2236      }
2237      if (n)
2238        coinFactorizationA_->emptyRows(n,which);
2239      delete [] which;
2240    }
2241  }
2242#endif
2243#ifdef CLP_FACTORIZATION_INSTRUMENT
2244  factorization_instrument(2);
2245#endif
2246  return coinFactorizationA_->status();
2247}
2248/* Replaces one Column in basis,
2249   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
2250   If checkBeforeModifying is true will do all accuracy checks
2251   before modifying factorization.  Whether to set this depends on
2252   speed considerations.  You could just do this on first iteration
2253   after factorization and thereafter re-factorize
2254   partial update already in U */
2255int 
2256ClpFactorization::replaceColumn ( const ClpSimplex * model, 
2257                                  CoinIndexedVector * regionSparse,
2258                                  CoinIndexedVector * tableauColumn,
2259                                  int pivotRow,
2260                                  double pivotCheck ,
2261                                  bool checkBeforeModifying,
2262                                       double acceptablePivot)
2263{
2264#ifndef SLIM_CLP
2265  if (!networkBasis_) {
2266#endif
2267#ifdef CLP_FACTORIZATION_INSTRUMENT
2268    factorization_instrument(-1);
2269#endif
2270    int returnCode;
2271    // see if FT
2272    if (!coinFactorizationA_||coinFactorizationA_->forrestTomlin()) {
2273      if (coinFactorizationA_) {
2274        returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2275                                                        pivotRow,
2276                                                        pivotCheck,
2277                                                        checkBeforeModifying,
2278                                                        acceptablePivot);
2279      } else {
2280        bool tab = coinFactorizationB_->wantsTableauColumn();
2281        int numberIterations=model->numberIterations();
2282        coinFactorizationB_->setUsefulInformation(&numberIterations,1);
2283        returnCode= 
2284          coinFactorizationB_->replaceColumn(tab?tableauColumn:regionSparse,
2285                                             pivotRow,
2286                                             pivotCheck,
2287                                             checkBeforeModifying,
2288                                             acceptablePivot);
2289#ifdef CLP_DEBUG
2290        // check basic
2291        int numberRows = coinFactorizationB_->numberRows();
2292        CoinIndexedVector region1(2*numberRows);
2293        CoinIndexedVector region2A(2*numberRows);
2294        CoinIndexedVector region2B(2*numberRows);
2295        int iPivot;
2296        double * arrayB = region2B.denseVector();
2297        int * pivotVariable = model->pivotVariable();
2298        int i;
2299        for (iPivot=0;iPivot<numberRows;iPivot++) {
2300          int iSequence = pivotVariable[iPivot];
2301          if (iPivot==pivotRow)
2302            iSequence = model->sequenceIn();
2303          model->unpack(&region2B,iSequence);
2304          coinFactorizationB_->updateColumn(&region1,&region2B);
2305          assert (fabs(arrayB[iPivot]-1.0)<1.0e-4);
2306          arrayB[iPivot]=0.0;
2307          for (i=0;i<numberRows;i++)
2308            assert (fabs(arrayB[i])<1.0e-4);
2309          region2B.clear();
2310        }
2311#endif
2312      }
2313    } else {
2314      returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn,
2315                                              pivotRow,pivotCheck); // Note array
2316    }
2317#ifdef CLP_FACTORIZATION_INSTRUMENT
2318      factorization_instrument(3);
2319#endif
2320      return returnCode;
2321
2322#ifndef SLIM_CLP
2323  } else {
2324    if (doCheck) {
2325      int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2326                                                pivotRow,
2327                                                pivotCheck,
2328                                                          checkBeforeModifying,
2329                                                          acceptablePivot);
2330      networkBasis_->replaceColumn(regionSparse,
2331                                   pivotRow);
2332      return returnCode;
2333    } else {
2334      // increase number of pivots
2335      coinFactorizationA_->setPivots(coinFactorizationA_->pivots()+1);
2336      return networkBasis_->replaceColumn(regionSparse,
2337                                   pivotRow);
2338    }
2339  }
2340#endif
2341}
2342
2343/* Updates one column (FTRAN) from region2
2344   number returned is negative if no room
2345   region1 starts as zero and is zero at end */
2346int 
2347ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
2348                                   CoinIndexedVector * regionSparse2)
2349{
2350#ifdef CLP_DEBUG
2351  regionSparse->checkClear();
2352#endif
2353  if (!numberRows())
2354      return 0;
2355#ifndef SLIM_CLP
2356  if (!networkBasis_) {
2357#endif
2358#ifdef CLP_FACTORIZATION_INSTRUMENT
2359    factorization_instrument(-1);
2360#endif
2361    int returnCode;
2362    if (coinFactorizationA_) {
2363      coinFactorizationA_->setCollectStatistics(true);
2364      returnCode= coinFactorizationA_->updateColumnFT(regionSparse,
2365                                                           regionSparse2);
2366      coinFactorizationA_->setCollectStatistics(false);
2367    } else {
2368      returnCode= coinFactorizationB_->updateColumnFT(regionSparse,
2369                                                         regionSparse2);
2370    }
2371#ifdef CLP_FACTORIZATION_INSTRUMENT
2372    factorization_instrument(4);
2373#endif
2374    return returnCode;
2375#ifndef SLIM_CLP
2376  } else {
2377#ifdef CHECK_NETWORK
2378    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2379    double * check = new double[coinFactorizationA_->numberRows()];
2380    int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2381                                                       regionSparse2);
2382    networkBasis_->updateColumn(regionSparse,save,-1);
2383    int i;
2384    double * array = regionSparse2->denseVector();
2385    int * indices = regionSparse2->getIndices();
2386    int n=regionSparse2->getNumElements();
2387    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2388    double * array2 = save->denseVector();
2389    int * indices2 = save->getIndices();
2390    int n2=save->getNumElements();
2391    assert (n==n2);
2392    if (save->packedMode()) {
2393      for (i=0;i<n;i++) {
2394        check[indices[i]]=array[i];
2395      }
2396      for (i=0;i<n;i++) {
2397        double value2 = array2[i];
2398        assert (check[indices2[i]]==value2);
2399      }
2400    } else {
2401      int numberRows = coinFactorizationA_->numberRows();
2402      for (i=0;i<numberRows;i++) {
2403        double value1 = array[i];
2404        double value2 = array2[i];
2405        assert (value1==value2);
2406      }
2407    }
2408    delete save;
2409    delete [] check;
2410    return returnCode;
2411#else
2412    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
2413    return 1;
2414#endif
2415  }
2416#endif
2417}
2418/* Updates one column (FTRAN) from region2
2419   number returned is negative if no room
2420   region1 starts as zero and is zero at end */
2421int 
2422ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
2423                                 CoinIndexedVector * regionSparse2,
2424                                 bool noPermute) const
2425{
2426#ifdef CLP_DEBUG
2427  if (!noPermute)
2428    regionSparse->checkClear();
2429#endif
2430  if (!numberRows())
2431    return 0;
2432#ifndef SLIM_CLP
2433  if (!networkBasis_) {
2434#endif
2435#ifdef CLP_FACTORIZATION_INSTRUMENT
2436    factorization_instrument(-1);
2437#endif
2438    int returnCode;
2439    if (coinFactorizationA_) {
2440      coinFactorizationA_->setCollectStatistics(true);
2441      returnCode= coinFactorizationA_->updateColumn(regionSparse,
2442                                                         regionSparse2,
2443                                                         noPermute);
2444      coinFactorizationA_->setCollectStatistics(false);
2445    } else {
2446      returnCode= coinFactorizationB_->updateColumn(regionSparse,
2447                                                         regionSparse2,
2448                                                         noPermute);
2449    }
2450#ifdef CLP_FACTORIZATION_INSTRUMENT
2451    factorization_instrument(5);
2452#endif
2453    //#define PRINT_VECTOR
2454#ifdef PRINT_VECTOR
2455    printf("Update\n");
2456    regionSparse2->print();
2457#endif
2458    return returnCode;
2459#ifndef SLIM_CLP
2460  } else {
2461#ifdef CHECK_NETWORK
2462    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2463    double * check = new double[coinFactorizationA_->numberRows()];
2464    int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2465                                                     regionSparse2,
2466                                                     noPermute);
2467    networkBasis_->updateColumn(regionSparse,save,-1);
2468    int i;
2469    double * array = regionSparse2->denseVector();
2470    int * indices = regionSparse2->getIndices();
2471    int n=regionSparse2->getNumElements();
2472    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2473    double * array2 = save->denseVector();
2474    int * indices2 = save->getIndices();
2475    int n2=save->getNumElements();
2476    assert (n==n2);
2477    if (save->packedMode()) {
2478      for (i=0;i<n;i++) {
2479        check[indices[i]]=array[i];
2480      }
2481      for (i=0;i<n;i++) {
2482        double value2 = array2[i];
2483        assert (check[indices2[i]]==value2);
2484      }
2485    } else {
2486      int numberRows = coinFactorizationA_->numberRows();
2487      for (i=0;i<numberRows;i++) {
2488        double value1 = array[i];
2489        double value2 = array2[i];
2490        assert (value1==value2);
2491      }
2492    }
2493    delete save;
2494    delete [] check;
2495    return returnCode;
2496#else
2497    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
2498    return 1;
2499#endif
2500  }
2501#endif
2502}
2503/* Updates one column (FTRAN) from region2
2504   Tries to do FT update
2505   number returned is negative if no room.
2506   Also updates region3
2507   region1 starts as zero and is zero at end */
2508int 
2509ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
2510                                       CoinIndexedVector * regionSparse2,
2511                                       CoinIndexedVector * regionSparse3,
2512                                       bool noPermuteRegion3)
2513{
2514#ifdef CLP_DEBUG
2515  regionSparse1->checkClear();
2516#endif
2517  if (!numberRows())
2518    return 0;
2519  int returnCode=0;
2520#ifndef SLIM_CLP
2521  if (!networkBasis_) {
2522#endif
2523#ifdef CLP_FACTORIZATION_INSTRUMENT
2524    factorization_instrument(-1);
2525#endif
2526    if (coinFactorizationA_) {
2527      coinFactorizationA_->setCollectStatistics(true);
2528      if (coinFactorizationA_->spaceForForrestTomlin()) {
2529        assert (regionSparse2->packedMode());
2530        assert (!regionSparse3->packedMode());
2531        returnCode= coinFactorizationA_->updateTwoColumnsFT(regionSparse1,
2532                                                            regionSparse2,
2533                                                            regionSparse3,
2534                                                            noPermuteRegion3);
2535      } else {
2536        returnCode= coinFactorizationA_->updateColumnFT(regionSparse1,
2537                                                        regionSparse2);
2538        coinFactorizationA_->updateColumn(regionSparse1,
2539                                          regionSparse3,
2540                                          noPermuteRegion3);
2541      }
2542      coinFactorizationA_->setCollectStatistics(false);
2543    } else {
2544#if 0
2545      CoinSimpFactorization * fact =
2546        dynamic_cast< CoinSimpFactorization*>(coinFactorizationB_);
2547      if (!fact) {
2548        returnCode= coinFactorizationB_->updateColumnFT(regionSparse1,
2549                                                        regionSparse2);
2550        coinFactorizationB_->updateColumn(regionSparse1,
2551                                          regionSparse3,
2552                                          noPermuteRegion3);
2553      } else {
2554        returnCode= fact->updateTwoColumnsFT(regionSparse1,
2555                                             regionSparse2,
2556                                             regionSparse3,
2557                                             noPermuteRegion3);
2558      }
2559#else
2560      returnCode= coinFactorizationB_->updateTwoColumnsFT(regionSparse1,
2561                                                          regionSparse2,
2562                                                          regionSparse3,
2563                                                          noPermuteRegion3);
2564#endif
2565    }
2566#ifdef CLP_FACTORIZATION_INSTRUMENT
2567    factorization_instrument(9);
2568#endif
2569#ifdef PRINT_VECTOR
2570    printf("UpdateTwoFT\n");
2571    regionSparse2->print();
2572    regionSparse3->print();
2573#endif
2574    return returnCode;
2575#ifndef SLIM_CLP
2576  } else {
2577    returnCode=updateColumnFT(regionSparse1,regionSparse2);
2578    updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
2579  }
2580#endif
2581  return returnCode;
2582}
2583/* Updates one column (FTRAN) from region2
2584   number returned is negative if no room
2585   region1 starts as zero and is zero at end */
2586int 
2587ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
2588                                 CoinIndexedVector * regionSparse2,
2589                                 bool noPermute) const
2590{
2591  if (!noPermute)
2592    regionSparse->checkClear();
2593  if (!coinFactorizationA_->numberRows())
2594    return 0;
2595  coinFactorizationA_->setCollectStatistics(false);
2596  int returnCode= coinFactorizationA_->updateColumn(regionSparse,
2597                                                   regionSparse2,
2598                                                   noPermute);
2599  return returnCode;
2600}
2601/* Updates one column (BTRAN) from region2
2602   region1 starts as zero and is zero at end */
2603int 
2604ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
2605                          CoinIndexedVector * regionSparse2) const
2606{
2607  if (!numberRows())
2608    return 0;
2609#ifndef SLIM_CLP
2610  if (!networkBasis_) {
2611#endif
2612#ifdef CLP_FACTORIZATION_INSTRUMENT
2613    factorization_instrument(-1);
2614#endif
2615    int returnCode;
2616
2617    if (coinFactorizationA_) {
2618      coinFactorizationA_->setCollectStatistics(true);
2619      returnCode =  coinFactorizationA_->updateColumnTranspose(regionSparse,
2620                                                    regionSparse2);
2621      coinFactorizationA_->setCollectStatistics(false);
2622    } else {
2623      returnCode= coinFactorizationB_->updateColumnTranspose(regionSparse,
2624                                                         regionSparse2);
2625    }
2626#ifdef CLP_FACTORIZATION_INSTRUMENT
2627    factorization_instrument(6);
2628#endif
2629#ifdef PRINT_VECTOR
2630    printf("UpdateTranspose\n");
2631    regionSparse2->print();
2632#endif
2633    return returnCode;
2634#ifndef SLIM_CLP
2635  } else {
2636#ifdef CHECK_NETWORK
2637    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2638    double * check = new double[coinFactorizationA_->numberRows()];
2639    int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
2640                                                              regionSparse2);
2641    networkBasis_->updateColumnTranspose(regionSparse,save);
2642    int i;
2643    double * array = regionSparse2->denseVector();
2644    int * indices = regionSparse2->getIndices();
2645    int n=regionSparse2->getNumElements();
2646    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2647    double * array2 = save->denseVector();
2648    int * indices2 = save->getIndices();
2649    int n2=save->getNumElements();
2650    assert (n==n2);
2651    if (save->packedMode()) {
2652      for (i=0;i<n;i++) {
2653        check[indices[i]]=array[i];
2654      }
2655      for (i=0;i<n;i++) {
2656        double value2 = array2[i];
2657        assert (check[indices2[i]]==value2);
2658      }
2659    } else {
2660      int numberRows = coinFactorizationA_->numberRows();
2661      for (i=0;i<numberRows;i++) {
2662        double value1 = array[i];
2663        double value2 = array2[i];
2664        assert (value1==value2);
2665      }
2666    }
2667    delete save;
2668    delete [] check;
2669    return returnCode;
2670#else
2671    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
2672#endif
2673  }
2674#endif
2675}
2676/* makes a row copy of L for speed and to allow very sparse problems */
2677void 
2678ClpFactorization::goSparse()
2679{
2680#ifndef SLIM_CLP
2681  if (!networkBasis_) {
2682#endif
2683    if (coinFactorizationA_) {
2684#ifdef CLP_FACTORIZATION_INSTRUMENT
2685      factorization_instrument(-1);
2686#endif
2687      coinFactorizationA_->goSparse();
2688#ifdef CLP_FACTORIZATION_INSTRUMENT
2689      factorization_instrument(7);
2690#endif
2691    }
2692  }
2693}
2694// Cleans up i.e. gets rid of network basis
2695void 
2696ClpFactorization::cleanUp()
2697{
2698#ifndef SLIM_CLP
2699  delete networkBasis_;
2700  networkBasis_=NULL;
2701#endif
2702  if (coinFactorizationA_)
2703  coinFactorizationA_->resetStatistics();
2704}
2705/// Says whether to redo pivot order
2706bool 
2707ClpFactorization::needToReorder() const
2708{
2709#ifdef CHECK_NETWORK
2710  return true;
2711#endif
2712#ifndef SLIM_CLP
2713  if (!networkBasis_)
2714#endif
2715    return true;
2716#ifndef SLIM_CLP
2717  else
2718    return false;
2719#endif
2720}
2721// Get weighted row list
2722void
2723ClpFactorization::getWeights(int * weights) const
2724{
2725#ifdef CLP_FACTORIZATION_INSTRUMENT
2726  factorization_instrument(-1);
2727#endif
2728#ifndef SLIM_CLP
2729  if (networkBasis_) {
2730    // Network - just unit
2731    int numberRows = coinFactorizationA_->numberRows();
2732    for (int i=0;i<numberRows;i++) 
2733      weights[i]=1;
2734    return;
2735  }
2736#endif
2737  int * numberInRow = coinFactorizationA_->numberInRow();
2738  int * numberInColumn = coinFactorizationA_->numberInColumn();
2739  int * permuteBack = coinFactorizationA_->pivotColumnBack();
2740  int * indexRowU = coinFactorizationA_->indexRowU();
2741  const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
2742  const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
2743  int numberRows = coinFactorizationA_->numberRows();
2744  if (!startRowL||!coinFactorizationA_->numberInRow()) {
2745    int * temp = new int[numberRows];
2746    memset(temp,0,numberRows*sizeof(int));
2747    int i;
2748    for (i=0;i<numberRows;i++) {
2749      // one for pivot
2750      temp[i]++;
2751      CoinBigIndex j;
2752      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
2753        int iRow=indexRowU[j];
2754        temp[iRow]++;
2755      }
2756    }
2757    CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
2758    int * indexRowL = coinFactorizationA_->indexRowL();
2759    int numberL = coinFactorizationA_->numberL();
2760    CoinBigIndex baseL = coinFactorizationA_->baseL();
2761    for (i=baseL;i<baseL+numberL;i++) {
2762      CoinBigIndex j;
2763      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
2764        int iRow = indexRowL[j];
2765        temp[iRow]++;
2766      }
2767    }
2768    for (i=0;i<numberRows;i++) {
2769      int number = temp[i];
2770      int iPermute = permuteBack[i];
2771      weights[iPermute]=number;
2772    }
2773    delete [] temp;
2774  } else {
2775    int i;
2776    for (i=0;i<numberRows;i++) {
2777      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
2778      //number = startRowL[i+1]-startRowL[i]+1;
2779      //number = numberInRow[i]+1;
2780      int iPermute = permuteBack[i];
2781      weights[iPermute]=number;
2782    }
2783  }
2784#ifdef CLP_FACTORIZATION_INSTRUMENT
2785  factorization_instrument(8);
2786#endif
2787}
2788// Set tolerances to safer of existing and given
2789void 
2790ClpFactorization::saferTolerances (  double zeroValue, 
2791                                    double pivotValue)
2792{
2793  // better to have small tolerance even if slower
2794  zeroTolerance(CoinMin(zeroTolerance(),zeroValue));
2795  // better to have large tolerance even if slower
2796  double newValue;
2797  if (pivotValue>0.0)
2798    newValue = pivotValue;
2799  else 
2800    newValue = -pivotTolerance()*pivotValue;
2801  pivotTolerance(CoinMin(CoinMax(pivotTolerance(),newValue),0.999));
2802}
2803// Sets factorization
2804void 
2805ClpFactorization::setFactorization(ClpFactorization & rhs)
2806{
2807  ClpFactorization::operator=(rhs);
2808}
2809#endif
Note: See TracBrowser for help on using the repository browser.