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

Last change on this file since 1287 was 1287, checked in by forrest, 13 years ago

changes for cbc event handler and multiple factorizations

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