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

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

changes for factorization and aux region

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