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

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

many changes to try and improve performance

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 76.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 CoinDenseFactorization();
1116  goDenseThreshold_ = -1;
1117}
1118
1119//-------------------------------------------------------------------
1120// Copy constructor
1121//-------------------------------------------------------------------
1122ClpFactorization::ClpFactorization (const ClpFactorization & rhs) 
1123{
1124#ifdef CLP_FACTORIZATION_INSTRUMENT
1125  factorization_instrument(-1);
1126#endif
1127#ifndef SLIM_CLP
1128  if (rhs.networkBasis_)
1129    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1130  else
1131    networkBasis_=NULL;
1132#endif
1133  if (rhs.coinFactorizationA_)
1134    coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
1135  else
1136    coinFactorizationA_=NULL;
1137  if (rhs.coinFactorizationB_)
1138    coinFactorizationB_ = new CoinDenseFactorization(*(rhs.coinFactorizationB_));
1139  else
1140    coinFactorizationB_=NULL;
1141  goDenseThreshold_ = rhs.goDenseThreshold_;
1142#ifdef CLP_FACTORIZATION_INSTRUMENT
1143  factorization_instrument(1);
1144#endif
1145}
1146
1147ClpFactorization::ClpFactorization (const CoinFactorization & rhs) 
1148{
1149#ifdef CLP_FACTORIZATION_INSTRUMENT
1150  factorization_instrument(-1);
1151#endif
1152#ifndef SLIM_CLP
1153  networkBasis_=NULL;
1154#endif
1155  coinFactorizationA_ = new CoinFactorization(rhs);
1156  coinFactorizationB_=NULL;
1157#ifdef CLP_FACTORIZATION_INSTRUMENT
1158  factorization_instrument(1);
1159  goDenseThreshold_ = -1;
1160#endif
1161}
1162
1163ClpFactorization::ClpFactorization (const CoinDenseFactorization & rhs) 
1164{
1165#ifdef CLP_FACTORIZATION_INSTRUMENT
1166  factorization_instrument(-1);
1167#endif
1168#ifndef SLIM_CLP
1169  networkBasis_=NULL;
1170#endif
1171  coinFactorizationA_ = NULL;
1172  coinFactorizationB_ = new CoinDenseFactorization(rhs);
1173  goDenseThreshold_ = -1;
1174#ifdef CLP_FACTORIZATION_INSTRUMENT
1175  factorization_instrument(1);
1176#endif
1177}
1178
1179//-------------------------------------------------------------------
1180// Destructor
1181//-------------------------------------------------------------------
1182ClpFactorization::~ClpFactorization () 
1183{
1184#ifndef SLIM_CLP
1185  delete networkBasis_;
1186#endif
1187  delete coinFactorizationA_;
1188  delete coinFactorizationB_;
1189}
1190
1191//----------------------------------------------------------------
1192// Assignment operator
1193//-------------------------------------------------------------------
1194ClpFactorization &
1195ClpFactorization::operator=(const ClpFactorization& rhs)
1196{
1197#ifdef CLP_FACTORIZATION_INSTRUMENT
1198  factorization_instrument(-1);
1199#endif
1200  if (this != &rhs) {
1201#ifndef SLIM_CLP
1202    delete networkBasis_;
1203    if (rhs.networkBasis_)
1204      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
1205    else
1206      networkBasis_=NULL;
1207#endif
1208    delete coinFactorizationA_;
1209    if (rhs.coinFactorizationA_)
1210      coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
1211    else
1212      coinFactorizationA_=NULL;
1213    delete coinFactorizationB_;
1214    if (rhs.coinFactorizationB_)
1215      coinFactorizationB_ = new CoinDenseFactorization(*(rhs.coinFactorizationB_));
1216    else
1217      coinFactorizationB_=NULL;
1218  }
1219#ifdef CLP_FACTORIZATION_INSTRUMENT
1220  factorization_instrument(1);
1221#endif
1222  return *this;
1223}
1224// Go over to dense code
1225void 
1226ClpFactorization::goDense() 
1227{
1228  delete coinFactorizationA_;
1229  delete coinFactorizationB_;
1230  coinFactorizationA_ = NULL;
1231  coinFactorizationB_ = new CoinDenseFactorization();
1232  //printf("going dense\n");
1233}
1234int 
1235ClpFactorization::factorize ( ClpSimplex * model,
1236                              int solveType, bool valuesPass)
1237{
1238  ClpMatrixBase * matrix = model->clpMatrix(); 
1239  int numberRows = model->numberRows();
1240  int numberColumns = model->numberColumns();
1241  if (!numberRows)
1242    return 0;
1243#ifdef CLP_FACTORIZATION_INSTRUMENT
1244  factorization_instrument(-1);
1245#endif
1246  if (coinFactorizationB_) {
1247    setStatus(-99);
1248    int * pivotVariable = model->pivotVariable();
1249    //returns 0 -okay, -1 singular, -2 too many in basis */
1250    while (status()<-98) {
1251     
1252      int i;
1253      int numberBasic=0;
1254      int numberRowBasic;
1255      // Move pivot variables across if they look good
1256      int * pivotTemp = model->rowArray(0)->getIndices();
1257      assert (!model->rowArray(0)->getNumElements());
1258      if (!matrix->rhsOffset(model)) {
1259        // Seems to prefer things in order so quickest
1260        // way is to go though like this
1261        for (i=0;i<numberRows;i++) {
1262          if (model->getRowStatus(i) == ClpSimplex::basic) 
1263            pivotTemp[numberBasic++]=i;
1264        }
1265        numberRowBasic=numberBasic;
1266        /* Put column basic variables into pivotVariable
1267           This is done by ClpMatrixBase to allow override for gub
1268        */
1269        matrix->generalExpanded(model,0,numberBasic);
1270      } else {
1271        // Long matrix - do a different way
1272        bool fullSearch=false;
1273        for (i=0;i<numberRows;i++) {
1274          int iPivot = pivotVariable[i];
1275          if (iPivot>=numberColumns) {
1276            pivotTemp[numberBasic++]=iPivot-numberColumns;
1277          }
1278        }
1279        numberRowBasic=numberBasic;
1280        for (i=0;i<numberRows;i++) {
1281          int iPivot = pivotVariable[i];
1282          if (iPivot<numberColumns) {
1283            if (iPivot>=0) {
1284              pivotTemp[numberBasic++]=iPivot;
1285            } else {
1286              // not full basis
1287              fullSearch=true;
1288              break;
1289            }
1290          }
1291        }
1292        if (fullSearch) {
1293          // do slow way
1294          numberBasic=0;
1295          for (i=0;i<numberRows;i++) {
1296            if (model->getRowStatus(i) == ClpSimplex::basic) 
1297              pivotTemp[numberBasic++]=i;
1298          }
1299          numberRowBasic=numberBasic;
1300          /* Put column basic variables into pivotVariable
1301             This is done by ClpMatrixBase to allow override for gub
1302          */
1303          matrix->generalExpanded(model,0,numberBasic);
1304        }
1305      }
1306      if (numberBasic>model->maximumBasic()) {
1307        // Take out some
1308        numberBasic=numberRowBasic;
1309        for (int i=0;i<numberColumns;i++) {
1310          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1311            if (numberBasic<numberRows)
1312              numberBasic++;
1313            else
1314              model->setColumnStatus(i,ClpSimplex::superBasic);
1315          }
1316        }
1317        numberBasic=numberRowBasic;
1318        matrix->generalExpanded(model,0,numberBasic);
1319      }
1320      CoinBigIndex numberElements=numberRowBasic;
1321     
1322      // compute how much in basis
1323      // can change for gub
1324      int numberColumnBasic = numberBasic-numberRowBasic;
1325     
1326      numberElements +=matrix->countBasis(model,
1327                                          pivotTemp+numberRowBasic, 
1328                                          numberRowBasic,
1329                                          numberColumnBasic);
1330      // Not needed for dense
1331      numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1332      coinFactorizationB_->getAreas ( numberRows,
1333                                      numberRowBasic+numberColumnBasic, numberElements,
1334                                      2 * numberElements );
1335      // Fill in counts so we can skip part of preProcess
1336      // This is NOT needed for dense but would be needed for later versions
1337      double * elementU;
1338      int * indexRowU;
1339      CoinBigIndex * startColumnU;
1340      int * numberInRow;
1341      int * numberInColumn;
1342      elementU = coinFactorizationB_->elements();
1343      indexRowU = coinFactorizationB_->indices();
1344      startColumnU = coinFactorizationB_->starts();
1345#ifndef COIN_FAST_CODE
1346      double slackValue;
1347      slackValue = coinFactorizationB_->slackValue();
1348#else
1349#define slackValue -1.0
1350#endif
1351      // At present we don't need counts
1352      numberInRow = coinFactorizationB_->intWorkArea();
1353      numberInColumn = numberInRow;
1354      CoinZeroN ( numberInRow, numberRows + 1 ); // just for valgrind
1355      for (i=0;i<numberRowBasic;i++) {
1356        int iRow = pivotTemp[i];
1357        // Change pivotTemp to correct sequence
1358        pivotTemp[i]=iRow+numberColumns;
1359        indexRowU[i]=iRow;
1360        startColumnU[i]=i;
1361        elementU[i]=slackValue;
1362      }
1363      startColumnU[numberRowBasic]=numberRowBasic;
1364      // can change for gub so redo
1365      numberColumnBasic = numberBasic-numberRowBasic;
1366      matrix->fillBasis(model, 
1367                        pivotTemp+numberRowBasic, 
1368                        numberColumnBasic,
1369                        indexRowU, 
1370                        startColumnU+numberRowBasic,
1371                        numberInRow,
1372                        numberInColumn+numberRowBasic,
1373                        elementU);
1374      // recompute number basic
1375      numberBasic = numberRowBasic+numberColumnBasic;
1376      for (i=numberBasic;i<numberRows;i++)
1377        pivotTemp[i]=-1; // mark not there
1378      if (numberBasic) 
1379        numberElements = startColumnU[numberBasic-1]
1380          +numberInColumn[numberBasic-1];
1381      else
1382        numberElements=0;
1383      coinFactorizationB_->preProcess ( );
1384      coinFactorizationB_->factor (  );
1385      // If we get here status is 0 or -1
1386      if (status() == 0&&numberBasic==numberRows) {
1387        coinFactorizationB_->postProcess(pivotTemp,pivotVariable);
1388      } else {
1389        // Change pivotTemp to be correct list
1390        coinFactorizationB_->makeNonSingular(pivotTemp,numberColumns);
1391        double * columnLower = model->lowerRegion();
1392        double * columnUpper = model->upperRegion();
1393        double * columnActivity = model->solutionRegion();
1394        double * rowLower = model->lowerRegion(0);
1395        double * rowUpper = model->upperRegion(0);
1396        double * rowActivity = model->solutionRegion(0);
1397        //redo basis - first take ALL out
1398        int iColumn;
1399        double largeValue = model->largeValue();
1400        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1401          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
1402            // take out
1403            if (!valuesPass) {
1404              double lower=columnLower[iColumn];
1405              double upper=columnUpper[iColumn];
1406              double value=columnActivity[iColumn];
1407              if (lower>-largeValue||upper<largeValue) {
1408                if (fabs(value-lower)<fabs(value-upper)) {
1409                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
1410                  columnActivity[iColumn]=lower;
1411                } else {
1412                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
1413                  columnActivity[iColumn]=upper;
1414                }
1415              } else {
1416                model->setColumnStatus(iColumn,ClpSimplex::isFree);
1417              }
1418            } else {
1419              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
1420            }
1421          }
1422        }
1423        int iRow;
1424        for (iRow=0;iRow<numberRows;iRow++) {
1425          if (model->getRowStatus(iRow)==ClpSimplex::basic) {
1426            // take out
1427            if (!valuesPass) {
1428              double lower=columnLower[iRow];
1429              double upper=columnUpper[iRow];
1430              double value=columnActivity[iRow];
1431              if (lower>-largeValue||upper<largeValue) {
1432                if (fabs(value-lower)<fabs(value-upper)) {
1433                  model->setRowStatus(iRow,ClpSimplex::atLowerBound);
1434                  columnActivity[iRow]=lower;
1435                } else {
1436                  model->setRowStatus(iRow,ClpSimplex::atUpperBound);
1437                  columnActivity[iRow]=upper;
1438                }
1439              } else {
1440                model->setRowStatus(iRow,ClpSimplex::isFree);
1441              }
1442            } else {
1443              model->setRowStatus(iRow,ClpSimplex::superBasic);
1444            }
1445          }
1446        }
1447        for (iRow=0;iRow<numberRows;iRow++) {
1448          int iSequence=pivotTemp[iRow];
1449          assert (iSequence>=0);
1450          // basic
1451          model->setColumnStatus(iSequence,ClpSimplex::basic);
1452        }
1453        // signal repeat
1454        coinFactorizationB_->setStatus(-99);
1455        // set fixed if they are
1456        for (iRow=0;iRow<numberRows;iRow++) {
1457          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
1458            if (rowLower[iRow]==rowUpper[iRow]) {
1459              rowActivity[iRow]=rowLower[iRow];
1460              model->setRowStatus(iRow,ClpSimplex::isFixed);
1461            }
1462          }
1463        }
1464        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1465          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
1466            if (columnLower[iColumn]==columnUpper[iColumn]) {
1467              columnActivity[iColumn]=columnLower[iColumn];
1468              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
1469            }
1470          }
1471        }
1472      } 
1473    }
1474#ifdef CLP_DEBUG
1475    // check basic
1476    CoinIndexedVector region1(2*numberRows);
1477    CoinIndexedVector region2B(2*numberRows);
1478    int iPivot;
1479    double * arrayB = region2B.denseVector();
1480    int i;
1481    for (iPivot=0;iPivot<numberRows;iPivot++) {
1482      int iSequence = pivotVariable[iPivot];
1483      model->unpack(&region2B,iSequence);
1484      coinFactorizationB_->updateColumn(&region1,&region2B);
1485      if (fabs(arrayB[iPivot]-1.0)<1.0e-4) {
1486        // OK?
1487        arrayB[iPivot]=0.0;
1488      } else {
1489        assert (fabs(arrayB[iPivot])<1.0e-4);
1490        for (i=0;i<numberRows;i++) {
1491          if (fabs(arrayB[i]-1.0)<1.0e-4)
1492            break;
1493        }
1494        assert (i<numberRows);
1495        printf("variable on row %d landed up on row %d\n",iPivot,i);
1496        arrayB[i]=0.0;
1497      }
1498      for (i=0;i<numberRows;i++)
1499        assert (fabs(arrayB[i])<1.0e-4);
1500      region2B.clear();
1501    }
1502#endif
1503#ifdef CLP_FACTORIZATION_INSTRUMENT
1504    factorization_instrument(2);
1505#endif
1506    return coinFactorizationB_->status();
1507  }
1508  // If too many compressions increase area
1509  if (coinFactorizationA_->pivots()>1&&coinFactorizationA_->numberCompressions()*10 > coinFactorizationA_->pivots()+10) {
1510    coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor()* 1.1);
1511  }
1512  //int numberPivots=coinFactorizationA_->pivots();
1513#if 0
1514  if (model->algorithm()>0)
1515    numberSave=-1;
1516#endif
1517#ifndef SLIM_CLP
1518  if (!networkBasis_||doCheck) {
1519#endif
1520    coinFactorizationA_->setStatus(-99);
1521    int * pivotVariable = model->pivotVariable();
1522    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
1523    while (coinFactorizationA_->status()<-98) {
1524     
1525      int i;
1526      int numberBasic=0;
1527      int numberRowBasic;
1528      // Move pivot variables across if they look good
1529      int * pivotTemp = model->rowArray(0)->getIndices();
1530      assert (!model->rowArray(0)->getNumElements());
1531      if (!matrix->rhsOffset(model)) {
1532#if 0
1533        if (numberSave>0) {
1534          int nStill=0;
1535          int nAtBound=0;
1536          int nZeroDual=0;
1537          CoinIndexedVector * array = model->rowArray(3);
1538          CoinIndexedVector * objArray = model->columnArray(1);
1539          array->clear();
1540          objArray->clear();
1541          double * cost = model->costRegion();
1542          double tolerance = model->primalTolerance();
1543          double offset=0.0;
1544          for (i=0;i<numberRows;i++) {
1545            int iPivot = pivotVariable[i];
1546            if (iPivot<numberColumns&&isDense(iPivot)) {
1547              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
1548                nStill++;
1549                double value=model->solutionRegion()[iPivot];
1550                double dual = model->dualRowSolution()[i];
1551                double lower=model->lowerRegion()[iPivot];
1552                double upper=model->upperRegion()[iPivot];
1553                ClpSimplex::Status status;
1554                if (fabs(value-lower)<tolerance) {
1555                  status=ClpSimplex::atLowerBound;
1556                  nAtBound++;
1557                } else if (fabs(value-upper)<tolerance) {
1558                  nAtBound++;
1559                  status=ClpSimplex::atUpperBound;
1560                } else if (value>lower&&value<upper) {
1561                  status=ClpSimplex::superBasic;
1562                } else {
1563                  status=ClpSimplex::basic;
1564                }
1565                if (status!=ClpSimplex::basic) {
1566                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
1567                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
1568                    model->setRowStatus(i,ClpSimplex::basic);
1569                    pivotVariable[i]=i+numberColumns;
1570                    model->dualRowSolution()[i]=0.0;
1571                    model->djRegion(0)[i]=0.0;
1572                    array->add(i,dual);
1573                    offset += dual*model->solutionRegion(0)[i];
1574                  }
1575                }
1576                if (fabs(dual)<1.0e-5)
1577                  nZeroDual++;
1578              }
1579            }
1580          }
1581          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
1582                 numberSave,nStill,nAtBound,nZeroDual,offset);
1583          if (array->getNumElements()) {
1584            // modify costs
1585            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
1586                                               objArray);
1587            array->clear();
1588            int n=objArray->getNumElements();
1589            int * indices = objArray->getIndices();
1590            double * elements = objArray->denseVector();
1591            for (i=0;i<n;i++) {
1592              int iColumn = indices[i];
1593              cost[iColumn] -= elements[iColumn];
1594              elements[iColumn]=0.0;
1595            }
1596            objArray->setNumElements(0);
1597          }
1598        }
1599#endif
1600        // Seems to prefer things in order so quickest
1601        // way is to go though like this
1602        for (i=0;i<numberRows;i++) {
1603          if (model->getRowStatus(i) == ClpSimplex::basic) 
1604            pivotTemp[numberBasic++]=i;
1605        }
1606        numberRowBasic=numberBasic;
1607        /* Put column basic variables into pivotVariable
1608           This is done by ClpMatrixBase to allow override for gub
1609        */
1610        matrix->generalExpanded(model,0,numberBasic);
1611      } else {
1612        // Long matrix - do a different way
1613        bool fullSearch=false;
1614        for (i=0;i<numberRows;i++) {
1615          int iPivot = pivotVariable[i];
1616          if (iPivot>=numberColumns) {
1617            pivotTemp[numberBasic++]=iPivot-numberColumns;
1618          }
1619        }
1620        numberRowBasic=numberBasic;
1621        for (i=0;i<numberRows;i++) {
1622          int iPivot = pivotVariable[i];
1623          if (iPivot<numberColumns) {
1624            if (iPivot>=0) {
1625              pivotTemp[numberBasic++]=iPivot;
1626            } else {
1627              // not full basis
1628              fullSearch=true;
1629              break;
1630            }
1631          }
1632        }
1633        if (fullSearch) {
1634          // do slow way
1635          numberBasic=0;
1636          for (i=0;i<numberRows;i++) {
1637            if (model->getRowStatus(i) == ClpSimplex::basic) 
1638              pivotTemp[numberBasic++]=i;
1639          }
1640          numberRowBasic=numberBasic;
1641          /* Put column basic variables into pivotVariable
1642             This is done by ClpMatrixBase to allow override for gub
1643          */
1644          matrix->generalExpanded(model,0,numberBasic);
1645        }
1646      }
1647      if (numberBasic>model->maximumBasic()) {
1648#if 0 // ndef NDEBUG
1649        printf("%d basic - should only be %d\n",
1650               numberBasic,numberRows);
1651#endif
1652        // Take out some
1653        numberBasic=numberRowBasic;
1654        for (int i=0;i<numberColumns;i++) {
1655          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1656            if (numberBasic<numberRows)
1657              numberBasic++;
1658            else
1659              model->setColumnStatus(i,ClpSimplex::superBasic);
1660          }
1661        }
1662        numberBasic=numberRowBasic;
1663        matrix->generalExpanded(model,0,numberBasic);
1664      }
1665#ifndef SLIM_CLP
1666      // see if matrix a network
1667#ifndef NO_RTTI
1668      ClpNetworkMatrix* networkMatrix =
1669        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
1670#else
1671      ClpNetworkMatrix* networkMatrix = NULL;
1672      if (model->clpMatrix()->type()==11)
1673        networkMatrix = 
1674        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
1675#endif
1676      // If network - still allow ordinary factorization first time for laziness
1677      if (networkMatrix)
1678        coinFactorizationA_->setBiasLU(0); // All to U if network
1679      //int saveMaximumPivots = maximumPivots();
1680      delete networkBasis_;
1681      networkBasis_ = NULL;
1682      if (networkMatrix&&!doCheck)
1683        maximumPivots(1);
1684#endif
1685      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
1686      while (coinFactorizationA_->status()==-99) {
1687        // maybe for speed will be better to leave as many regions as possible
1688        coinFactorizationA_->gutsOfDestructor();
1689        coinFactorizationA_->gutsOfInitialize(2);
1690        CoinBigIndex numberElements=numberRowBasic;
1691
1692        // compute how much in basis
1693
1694        int i;
1695        // can change for gub
1696        int numberColumnBasic = numberBasic-numberRowBasic;
1697
1698        numberElements +=matrix->countBasis(model,
1699                                           pivotTemp+numberRowBasic, 
1700                                           numberRowBasic,
1701                                            numberColumnBasic);
1702        // and recompute as network side say different
1703        if (model->numberIterations())
1704          numberRowBasic = numberBasic - numberColumnBasic;
1705        numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1706        coinFactorizationA_->getAreas ( numberRows,
1707                   numberRowBasic+numberColumnBasic, numberElements,
1708                   2 * numberElements );
1709        //fill
1710        // Fill in counts so we can skip part of preProcess
1711        int * numberInRow = coinFactorizationA_->numberInRow();
1712        int * numberInColumn = coinFactorizationA_->numberInColumn();
1713        CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
1714        CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
1715        double * elementU = coinFactorizationA_->elementU();
1716        int * indexRowU = coinFactorizationA_->indexRowU();
1717        CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
1718#ifndef COIN_FAST_CODE
1719        double slackValue = coinFactorizationA_->slackValue();
1720#endif
1721        for (i=0;i<numberRowBasic;i++) {
1722          int iRow = pivotTemp[i];
1723          indexRowU[i]=iRow;
1724          startColumnU[i]=i;
1725          elementU[i]=slackValue;
1726          numberInRow[iRow]=1;
1727          numberInColumn[i]=1;
1728        }
1729        startColumnU[numberRowBasic]=numberRowBasic;
1730        // can change for gub so redo
1731        numberColumnBasic = numberBasic-numberRowBasic;
1732        matrix->fillBasis(model, 
1733                          pivotTemp+numberRowBasic, 
1734                          numberColumnBasic,
1735                          indexRowU, 
1736                          startColumnU+numberRowBasic,
1737                          numberInRow,
1738                          numberInColumn+numberRowBasic,
1739                          elementU);
1740#if 0
1741        {
1742          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
1743          for (int i=0;i<numberElements;i++)
1744            printf("row %d col %d value %g\n",indexRowU[i],indexColumnU_[i],
1745                   elementU[i]);
1746        }
1747#endif
1748        // recompute number basic
1749        numberBasic = numberRowBasic+numberColumnBasic;
1750        if (numberBasic) 
1751          numberElements = startColumnU[numberBasic-1]
1752            +numberInColumn[numberBasic-1];
1753        else
1754          numberElements=0;
1755        coinFactorizationA_->setNumberElementsU(numberElements);
1756        //saveFactorization("dump.d");
1757        if (coinFactorizationA_->biasLU()>=3||coinFactorizationA_->numberRows()!=coinFactorizationA_->numberColumns())
1758          coinFactorizationA_->preProcess ( 2 );
1759        else
1760          coinFactorizationA_->preProcess ( 3 ); // no row copy
1761        coinFactorizationA_->factor (  );
1762        if (coinFactorizationA_->status()==-99) {
1763          // get more memory
1764          coinFactorizationA_->areaFactor(2.0*coinFactorizationA_->areaFactor());
1765        } else if (coinFactorizationA_->status()==-1&&model->numberIterations()==0&&
1766                   coinFactorizationA_->denseThreshold()) {
1767          // Round again without dense
1768          coinFactorizationA_->setDenseThreshold(0);
1769          coinFactorizationA_->setStatus(-99);
1770        }
1771      }
1772      // If we get here status is 0 or -1
1773      if (coinFactorizationA_->status() == 0) {
1774        // We may need to tamper with order and redo - e.g. network with side
1775        int useNumberRows = numberRows;
1776        // **** we will also need to add test in dual steepest to do
1777        // as we do for network
1778        matrix->generalExpanded(model,12,useNumberRows);
1779        const int * permuteBack = coinFactorizationA_->permuteBack();
1780        const int * back = coinFactorizationA_->pivotColumnBack();
1781        //int * pivotTemp = pivotColumn_.array();
1782        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
1783        // Redo pivot order
1784        for (i=0;i<numberRowBasic;i++) {
1785          int k = pivotTemp[i];
1786          // so rowIsBasic[k] would be permuteBack[back[i]]
1787          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
1788        }
1789        for (;i<useNumberRows;i++) {
1790          int k = pivotTemp[i];
1791          // so rowIsBasic[k] would be permuteBack[back[i]]
1792          pivotVariable[permuteBack[back[i]]]=k;
1793        }
1794#if 0
1795        if (numberSave>=0) {
1796          numberSave=numberDense_;
1797          memset(saveList,0,((coinFactorizationA_->numberRows()+31)>>5)*sizeof(int));
1798          for (i=coinFactorizationA_->numberRows()-numberSave;i<coinFactorizationA_->numberRows();i++) {
1799            int k=pivotTemp[pivotColumn_.array()[i]];
1800            setDense(k);
1801          }
1802        }
1803#endif
1804        // Set up permutation vector
1805        // these arrays start off as copies of permute
1806        // (and we could use permute_ instead of pivotColumn (not back though))
1807        ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn()  );
1808        ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack()  );
1809#ifndef SLIM_CLP
1810        if (networkMatrix) {
1811          maximumPivots(CoinMax(2000,maximumPivots()));
1812          // redo arrays
1813          for (int iRow=0;iRow<4;iRow++) {
1814            int length =model->numberRows()+maximumPivots();
1815            if (iRow==3||model->objectiveAsObject()->type()>1)
1816              length += model->numberColumns();
1817            model->rowArray(iRow)->reserve(length);
1818          }
1819          // create network factorization
1820          if (doCheck)
1821            delete networkBasis_; // temp
1822          networkBasis_ = new ClpNetworkBasis(model,coinFactorizationA_->numberRows(),
1823                                              coinFactorizationA_->pivotRegion(),
1824                                              coinFactorizationA_->permuteBack(),
1825                                              coinFactorizationA_->startColumnU(),
1826                                              coinFactorizationA_->numberInColumn(),
1827                                              coinFactorizationA_->indexRowU(),
1828                                              coinFactorizationA_->elementU());
1829          // kill off arrays in ordinary factorization
1830          if (!doCheck) {
1831            coinFactorizationA_->gutsOfDestructor();
1832            // but make sure coinFactorizationA_->numberRows() set
1833            coinFactorizationA_->setNumberRows(model->numberRows());
1834            coinFactorizationA_->setStatus(0);
1835#if 0
1836            // but put back permute arrays so odd things will work
1837            int numberRows = model->numberRows();
1838            pivotColumnBack_ = new int [numberRows];
1839            permute_ = new int [numberRows];
1840            int i;
1841            for (i=0;i<numberRows;i++) {
1842              pivotColumnBack_[i]=i;
1843              permute_[i]=i;
1844            }
1845#endif
1846          }
1847        } else {
1848#endif
1849          // See if worth going sparse and when
1850          coinFactorizationA_->checkSparse();
1851#ifndef SLIM_CLP
1852        }
1853#endif
1854      } else if (coinFactorizationA_->status() == -1&&(solveType==0||solveType==2)) {
1855        // This needs redoing as it was merged coding - does not need array
1856        int numberTotal = numberRows+numberColumns;
1857        int * isBasic = new int [numberTotal];
1858        int * rowIsBasic = isBasic+numberColumns;
1859        int * columnIsBasic = isBasic;
1860        for (i=0;i<numberTotal;i++) 
1861          isBasic[i]=-1;
1862        for (i=0;i<numberRowBasic;i++) {
1863          int iRow = pivotTemp[i];
1864          rowIsBasic[iRow]=1;
1865        }
1866        for (;i<numberBasic;i++) {
1867          int iColumn = pivotTemp[i];
1868          columnIsBasic[iColumn]=1;
1869        }
1870        numberBasic=0;
1871        for (i=0;i<numberRows;i++) 
1872          pivotVariable[i]=-1;
1873        // mark as basic or non basic
1874        const int * pivotColumn = coinFactorizationA_->pivotColumn();
1875        for (i=0;i<numberRows;i++) {
1876          if (rowIsBasic[i]>=0) {
1877            if (pivotColumn[numberBasic]>=0) {
1878              rowIsBasic[i]=pivotColumn[numberBasic];
1879            } else {
1880              rowIsBasic[i]=-1;
1881              model->setRowStatus(i,ClpSimplex::superBasic);
1882            }
1883            numberBasic++;
1884          }
1885        }
1886        for (i=0;i<numberColumns;i++) {
1887          if (columnIsBasic[i]>=0) {
1888            if (pivotColumn[numberBasic]>=0) 
1889              columnIsBasic[i]=pivotColumn[numberBasic];
1890            else
1891              columnIsBasic[i]=-1;
1892            numberBasic++;
1893          }
1894        }
1895        // leave pivotVariable in useful form for cleaning basis
1896        int * pivotVariable = model->pivotVariable();
1897        for (i=0;i<numberRows;i++) {
1898          pivotVariable[i]=-1;
1899        }
1900
1901        for (i=0;i<numberRows;i++) {
1902          if (model->getRowStatus(i) == ClpSimplex::basic) {
1903            int iPivot = rowIsBasic[i];
1904            if (iPivot>=0) 
1905              pivotVariable[iPivot]=i+numberColumns;
1906          }
1907        }
1908        for (i=0;i<numberColumns;i++) {
1909          if (model->getColumnStatus(i) == ClpSimplex::basic) {
1910            int iPivot = columnIsBasic[i];
1911            if (iPivot>=0) 
1912              pivotVariable[iPivot]=i;
1913          }
1914        }
1915        delete [] isBasic;
1916        double * columnLower = model->lowerRegion();
1917        double * columnUpper = model->upperRegion();
1918        double * columnActivity = model->solutionRegion();
1919        double * rowLower = model->lowerRegion(0);
1920        double * rowUpper = model->upperRegion(0);
1921        double * rowActivity = model->solutionRegion(0);
1922        //redo basis - first take ALL columns out
1923        int iColumn;
1924        double largeValue = model->largeValue();
1925        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1926          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
1927            // take out
1928            if (!valuesPass) {
1929              double lower=columnLower[iColumn];
1930              double upper=columnUpper[iColumn];
1931              double value=columnActivity[iColumn];
1932              if (lower>-largeValue||upper<largeValue) {
1933                if (fabs(value-lower)<fabs(value-upper)) {
1934                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
1935                  columnActivity[iColumn]=lower;
1936                } else {
1937                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
1938                  columnActivity[iColumn]=upper;
1939                }
1940              } else {
1941                model->setColumnStatus(iColumn,ClpSimplex::isFree);
1942              }
1943            } else {
1944              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
1945            }
1946          }
1947        }
1948        int iRow;
1949        for (iRow=0;iRow<numberRows;iRow++) {
1950          int iSequence=pivotVariable[iRow];
1951          if (iSequence>=0) {
1952            // basic
1953            if (iSequence>=numberColumns) {
1954              // slack in - leave
1955              //assert (iSequence-numberColumns==iRow);
1956            } else {
1957              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
1958              // put back structural
1959              model->setColumnStatus(iSequence,ClpSimplex::basic);
1960            }
1961          } else {
1962            // put in slack
1963            model->setRowStatus(iRow,ClpSimplex::basic);
1964          }
1965        }
1966        // Put back any key variables for gub
1967        int dummy;
1968        matrix->generalExpanded(model,1,dummy);
1969        // signal repeat
1970        coinFactorizationA_->setStatus(-99);
1971        // set fixed if they are
1972        for (iRow=0;iRow<numberRows;iRow++) {
1973          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
1974            if (rowLower[iRow]==rowUpper[iRow]) {
1975              rowActivity[iRow]=rowLower[iRow];
1976              model->setRowStatus(iRow,ClpSimplex::isFixed);
1977            }
1978          }
1979        }
1980        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1981          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
1982            if (columnLower[iColumn]==columnUpper[iColumn]) {
1983              columnActivity[iColumn]=columnLower[iColumn];
1984              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
1985            }
1986          }
1987        }
1988      } 
1989    }
1990#ifndef SLIM_CLP
1991  } else {
1992    // network - fake factorization - do nothing
1993    coinFactorizationA_->setStatus(0);
1994    coinFactorizationA_->setPivots(0);
1995  }
1996#endif
1997#ifndef SLIM_CLP
1998  if (!coinFactorizationA_->status()) {
1999    // take out part if quadratic
2000    if (model->algorithm()==2) {
2001      ClpObjective * obj = model->objectiveAsObject();
2002      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2003      assert (quadraticObj);
2004      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2005      int numberXColumns = quadratic->getNumCols();
2006      assert (numberXColumns<numberColumns);
2007      int base = numberColumns-numberXColumns;
2008      int * which = new int [numberXColumns];
2009      int * pivotVariable = model->pivotVariable();
2010      int * permute = pivotColumn();
2011      int i;
2012      int n=0;
2013      for (i=0;i<numberRows;i++) {
2014        int iSj = pivotVariable[i]-base;
2015        if (iSj>=0&&iSj<numberXColumns) 
2016          which[n++]=permute[i];
2017      }
2018      if (n)
2019        coinFactorizationA_->emptyRows(n,which);
2020      delete [] which;
2021    }
2022  }
2023#endif
2024#ifdef CLP_FACTORIZATION_INSTRUMENT
2025  factorization_instrument(2);
2026#endif
2027  return coinFactorizationA_->status();
2028}
2029/* Replaces one Column to basis,
2030   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
2031   If checkBeforeModifying is true will do all accuracy checks
2032   before modifying factorization.  Whether to set this depends on
2033   speed considerations.  You could just do this on first iteration
2034   after factorization and thereafter re-factorize
2035   partial update already in U */
2036int 
2037ClpFactorization::replaceColumn ( const ClpSimplex * model, 
2038                                  CoinIndexedVector * regionSparse,
2039                                  CoinIndexedVector * tableauColumn,
2040                                  int pivotRow,
2041                                  double pivotCheck ,
2042                                  bool checkBeforeModifying)
2043{
2044#ifndef SLIM_CLP
2045  if (!networkBasis_) {
2046#endif
2047#ifdef CLP_FACTORIZATION_INSTRUMENT
2048    factorization_instrument(-1);
2049#endif
2050    int returnCode;
2051    // see if FT
2052    if (!coinFactorizationA_||coinFactorizationA_->forrestTomlin()) {
2053      if (coinFactorizationA_) {
2054        returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2055                                                        pivotRow,
2056                                                        pivotCheck,
2057                                                        checkBeforeModifying);
2058      } else {
2059        returnCode= coinFactorizationB_->replaceColumn(tableauColumn,
2060                                              pivotRow,
2061                                              pivotCheck,
2062                                              checkBeforeModifying);
2063#ifdef CLP_DEBUG
2064        // check basic
2065        int numberRows = coinFactorizationB_->numberRows();
2066        CoinIndexedVector region1(2*numberRows);
2067        CoinIndexedVector region2A(2*numberRows);
2068        CoinIndexedVector region2B(2*numberRows);
2069        int iPivot;
2070        double * arrayB = region2B.denseVector();
2071        int * pivotVariable = model->pivotVariable();
2072        int i;
2073        for (iPivot=0;iPivot<numberRows;iPivot++) {
2074          int iSequence = pivotVariable[iPivot];
2075          if (iPivot==pivotRow)
2076            iSequence = model->sequenceIn();
2077          model->unpack(&region2B,iSequence);
2078          coinFactorizationB_->updateColumn(&region1,&region2B);
2079          assert (fabs(arrayB[iPivot]-1.0)<1.0e-4);
2080          arrayB[iPivot]=0.0;
2081          for (i=0;i<numberRows;i++)
2082            assert (fabs(arrayB[i])<1.0e-4);
2083          region2B.clear();
2084        }
2085#endif
2086      }
2087    } else {
2088      returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn,
2089                                              pivotRow,pivotCheck); // Note array
2090    }
2091#ifdef CLP_FACTORIZATION_INSTRUMENT
2092      factorization_instrument(3);
2093#endif
2094      return returnCode;
2095
2096#ifndef SLIM_CLP
2097  } else {
2098    if (doCheck) {
2099      int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2100                                                pivotRow,
2101                                                pivotCheck,
2102                                                checkBeforeModifying);
2103      networkBasis_->replaceColumn(regionSparse,
2104                                   pivotRow);
2105      return returnCode;
2106    } else {
2107      // increase number of pivots
2108      coinFactorizationA_->setPivots(coinFactorizationA_->pivots()+1);
2109      return networkBasis_->replaceColumn(regionSparse,
2110                                   pivotRow);
2111    }
2112  }
2113#endif
2114}
2115
2116/* Updates one column (FTRAN) from region2
2117   number returned is negative if no room
2118   region1 starts as zero and is zero at end */
2119int 
2120ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
2121                                   CoinIndexedVector * regionSparse2)
2122{
2123#ifdef CLP_DEBUG
2124  regionSparse->checkClear();
2125#endif
2126  if (!numberRows())
2127      return 0;
2128#ifndef SLIM_CLP
2129  if (!networkBasis_) {
2130#endif
2131#ifdef CLP_FACTORIZATION_INSTRUMENT
2132    factorization_instrument(-1);
2133#endif
2134    int returnCode;
2135    if (coinFactorizationA_) {
2136      coinFactorizationA_->setCollectStatistics(true);
2137      returnCode= coinFactorizationA_->updateColumnFT(regionSparse,
2138                                                           regionSparse2);
2139      coinFactorizationA_->setCollectStatistics(false);
2140    } else {
2141      returnCode= coinFactorizationB_->updateColumn(regionSparse,
2142                                                         regionSparse2);
2143    }
2144#ifdef CLP_FACTORIZATION_INSTRUMENT
2145    factorization_instrument(4);
2146#endif
2147    return returnCode;
2148#ifndef SLIM_CLP
2149  } else {
2150#ifdef CHECK_NETWORK
2151    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2152    double * check = new double[coinFactorizationA_->numberRows()];
2153    int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2154                                                       regionSparse2);
2155    networkBasis_->updateColumn(regionSparse,save,-1);
2156    int i;
2157    double * array = regionSparse2->denseVector();
2158    int * indices = regionSparse2->getIndices();
2159    int n=regionSparse2->getNumElements();
2160    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2161    double * array2 = save->denseVector();
2162    int * indices2 = save->getIndices();
2163    int n2=save->getNumElements();
2164    assert (n==n2);
2165    if (save->packedMode()) {
2166      for (i=0;i<n;i++) {
2167        check[indices[i]]=array[i];
2168      }
2169      for (i=0;i<n;i++) {
2170        double value2 = array2[i];
2171        assert (check[indices2[i]]==value2);
2172      }
2173    } else {
2174      int numberRows = coinFactorizationA_->numberRows();
2175      for (i=0;i<numberRows;i++) {
2176        double value1 = array[i];
2177        double value2 = array2[i];
2178        assert (value1==value2);
2179      }
2180    }
2181    delete save;
2182    delete [] check;
2183    return returnCode;
2184#else
2185    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
2186    return 1;
2187#endif
2188  }
2189#endif
2190}
2191/* Updates one column (FTRAN) from region2
2192   number returned is negative if no room
2193   region1 starts as zero and is zero at end */
2194int 
2195ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
2196                                 CoinIndexedVector * regionSparse2,
2197                                 bool noPermute) const
2198{
2199#ifdef CLP_DEBUG
2200  if (!noPermute)
2201    regionSparse->checkClear();
2202#endif
2203  if (!numberRows())
2204    return 0;
2205#ifndef SLIM_CLP
2206  if (!networkBasis_) {
2207#endif
2208#ifdef CLP_FACTORIZATION_INSTRUMENT
2209    factorization_instrument(-1);
2210#endif
2211    int returnCode;
2212    if (coinFactorizationA_) {
2213      coinFactorizationA_->setCollectStatistics(true);
2214      returnCode= coinFactorizationA_->updateColumn(regionSparse,
2215                                                         regionSparse2,
2216                                                         noPermute);
2217      coinFactorizationA_->setCollectStatistics(false);
2218    } else {
2219      returnCode= coinFactorizationB_->updateColumn(regionSparse,
2220                                                         regionSparse2,
2221                                                         noPermute);
2222    }
2223#ifdef CLP_FACTORIZATION_INSTRUMENT
2224    factorization_instrument(5);
2225#endif
2226    return returnCode;
2227#ifndef SLIM_CLP
2228  } else {
2229#ifdef CHECK_NETWORK
2230    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2231    double * check = new double[coinFactorizationA_->numberRows()];
2232    int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2233                                                     regionSparse2,
2234                                                     noPermute);
2235    networkBasis_->updateColumn(regionSparse,save,-1);
2236    int i;
2237    double * array = regionSparse2->denseVector();
2238    int * indices = regionSparse2->getIndices();
2239    int n=regionSparse2->getNumElements();
2240    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2241    double * array2 = save->denseVector();
2242    int * indices2 = save->getIndices();
2243    int n2=save->getNumElements();
2244    assert (n==n2);
2245    if (save->packedMode()) {
2246      for (i=0;i<n;i++) {
2247        check[indices[i]]=array[i];
2248      }
2249      for (i=0;i<n;i++) {
2250        double value2 = array2[i];
2251        assert (check[indices2[i]]==value2);
2252      }
2253    } else {
2254      int numberRows = coinFactorizationA_->numberRows();
2255      for (i=0;i<numberRows;i++) {
2256        double value1 = array[i];
2257        double value2 = array2[i];
2258        assert (value1==value2);
2259      }
2260    }
2261    delete save;
2262    delete [] check;
2263    return returnCode;
2264#else
2265    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
2266    return 1;
2267#endif
2268  }
2269#endif
2270}
2271/* Updates one column (FTRAN) from region2
2272   Tries to do FT update
2273   number returned is negative if no room.
2274   Also updates region3
2275   region1 starts as zero and is zero at end */
2276int 
2277ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
2278                                       CoinIndexedVector * regionSparse2,
2279                                       CoinIndexedVector * regionSparse3,
2280                                       bool noPermuteRegion3)
2281{
2282#ifdef CLP_DEBUG
2283  regionSparse1->checkClear();
2284#endif
2285  if (!numberRows())
2286    return 0;
2287  int returnCode=0;
2288#ifndef SLIM_CLP
2289  if (!networkBasis_) {
2290#endif
2291#ifdef CLP_FACTORIZATION_INSTRUMENT
2292    factorization_instrument(-1);
2293#endif
2294    if (coinFactorizationA_) {
2295      coinFactorizationA_->setCollectStatistics(true);
2296      if (coinFactorizationA_->spaceForForrestTomlin()) {
2297        assert (regionSparse2->packedMode());
2298        assert (!regionSparse3->packedMode());
2299        returnCode= coinFactorizationA_->updateTwoColumnsFT(regionSparse1,
2300                                                            regionSparse2,
2301                                                            regionSparse3,
2302                                                            noPermuteRegion3);
2303      } else {
2304        returnCode= coinFactorizationA_->updateColumnFT(regionSparse1,
2305                                                        regionSparse2);
2306        coinFactorizationA_->updateColumn(regionSparse1,
2307                                          regionSparse3,
2308                                          noPermuteRegion3);
2309      }
2310      coinFactorizationA_->setCollectStatistics(false);
2311    } else {
2312      returnCode= coinFactorizationB_->updateColumn(regionSparse1,
2313                                                     regionSparse2);
2314      coinFactorizationB_->updateColumn(regionSparse1,
2315                                        regionSparse3,
2316                                        noPermuteRegion3);
2317    }
2318#ifdef CLP_FACTORIZATION_INSTRUMENT
2319    factorization_instrument(9);
2320#endif
2321    return returnCode;
2322#ifndef SLIM_CLP
2323  } else {
2324    returnCode=updateColumnFT(regionSparse1,regionSparse2);
2325    updateColumn(regionSparse1,regionSparse3,noPermuteRegion3);
2326  }
2327#endif
2328  return returnCode;
2329}
2330/* Updates one column (FTRAN) from region2
2331   number returned is negative if no room
2332   region1 starts as zero and is zero at end */
2333int 
2334ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
2335                                 CoinIndexedVector * regionSparse2,
2336                                 bool noPermute) const
2337{
2338  if (!noPermute)
2339    regionSparse->checkClear();
2340  if (!coinFactorizationA_->numberRows())
2341    return 0;
2342  coinFactorizationA_->setCollectStatistics(false);
2343  int returnCode= coinFactorizationA_->updateColumn(regionSparse,
2344                                                   regionSparse2,
2345                                                   noPermute);
2346  return returnCode;
2347}
2348/* Updates one column (BTRAN) from region2
2349   region1 starts as zero and is zero at end */
2350int 
2351ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
2352                          CoinIndexedVector * regionSparse2) const
2353{
2354  if (!numberRows())
2355    return 0;
2356#ifndef SLIM_CLP
2357  if (!networkBasis_) {
2358#endif
2359#ifdef CLP_FACTORIZATION_INSTRUMENT
2360    factorization_instrument(-1);
2361#endif
2362    int returnCode;
2363    if (coinFactorizationA_) {
2364      coinFactorizationA_->setCollectStatistics(true);
2365      returnCode =  coinFactorizationA_->updateColumnTranspose(regionSparse,
2366                                                    regionSparse2);
2367      coinFactorizationA_->setCollectStatistics(false);
2368    } else {
2369      returnCode= coinFactorizationB_->updateColumnTranspose(regionSparse,
2370                                                         regionSparse2);
2371    }
2372#ifdef CLP_FACTORIZATION_INSTRUMENT
2373    factorization_instrument(6);
2374#endif
2375    return returnCode;
2376#ifndef SLIM_CLP
2377  } else {
2378#ifdef CHECK_NETWORK
2379    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2380    double * check = new double[coinFactorizationA_->numberRows()];
2381    int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
2382                                                              regionSparse2);
2383    networkBasis_->updateColumnTranspose(regionSparse,save);
2384    int i;
2385    double * array = regionSparse2->denseVector();
2386    int * indices = regionSparse2->getIndices();
2387    int n=regionSparse2->getNumElements();
2388    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
2389    double * array2 = save->denseVector();
2390    int * indices2 = save->getIndices();
2391    int n2=save->getNumElements();
2392    assert (n==n2);
2393    if (save->packedMode()) {
2394      for (i=0;i<n;i++) {
2395        check[indices[i]]=array[i];
2396      }
2397      for (i=0;i<n;i++) {
2398        double value2 = array2[i];
2399        assert (check[indices2[i]]==value2);
2400      }
2401    } else {
2402      int numberRows = coinFactorizationA_->numberRows();
2403      for (i=0;i<numberRows;i++) {
2404        double value1 = array[i];
2405        double value2 = array2[i];
2406        assert (value1==value2);
2407      }
2408    }
2409    delete save;
2410    delete [] check;
2411    return returnCode;
2412#else
2413    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
2414#endif
2415  }
2416#endif
2417}
2418/* makes a row copy of L for speed and to allow very sparse problems */
2419void 
2420ClpFactorization::goSparse()
2421{
2422#ifndef SLIM_CLP
2423  if (!networkBasis_) {
2424#endif
2425    if (coinFactorizationA_) {
2426#ifdef CLP_FACTORIZATION_INSTRUMENT
2427      factorization_instrument(-1);
2428#endif
2429      coinFactorizationA_->goSparse();
2430#ifdef CLP_FACTORIZATION_INSTRUMENT
2431      factorization_instrument(7);
2432#endif
2433    }
2434  }
2435}
2436// Cleans up i.e. gets rid of network basis
2437void 
2438ClpFactorization::cleanUp()
2439{
2440#ifndef SLIM_CLP
2441  delete networkBasis_;
2442  networkBasis_=NULL;
2443#endif
2444  if (coinFactorizationA_)
2445  coinFactorizationA_->resetStatistics();
2446}
2447/// Says whether to redo pivot order
2448bool 
2449ClpFactorization::needToReorder() const
2450{
2451#ifdef CHECK_NETWORK
2452  return true;
2453#endif
2454#ifndef SLIM_CLP
2455  if (!networkBasis_)
2456#endif
2457    return true;
2458#ifndef SLIM_CLP
2459  else
2460    return false;
2461#endif
2462}
2463// Get weighted row list
2464void
2465ClpFactorization::getWeights(int * weights) const
2466{
2467#ifdef CLP_FACTORIZATION_INSTRUMENT
2468  factorization_instrument(-1);
2469#endif
2470#ifndef SLIM_CLP
2471  if (networkBasis_) {
2472    // Network - just unit
2473    int numberRows = coinFactorizationA_->numberRows();
2474    for (int i=0;i<numberRows;i++) 
2475      weights[i]=1;
2476    return;
2477  }
2478#endif
2479  int * numberInRow = coinFactorizationA_->numberInRow();
2480  int * numberInColumn = coinFactorizationA_->numberInColumn();
2481  int * permuteBack = coinFactorizationA_->pivotColumnBack();
2482  int * indexRowU = coinFactorizationA_->indexRowU();
2483  const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
2484  const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
2485  int numberRows = coinFactorizationA_->numberRows();
2486  if (!startRowL||!coinFactorizationA_->numberInRow()) {
2487    int * temp = new int[numberRows];
2488    memset(temp,0,numberRows*sizeof(int));
2489    int i;
2490    for (i=0;i<numberRows;i++) {
2491      // one for pivot
2492      temp[i]++;
2493      CoinBigIndex j;
2494      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
2495        int iRow=indexRowU[j];
2496        temp[iRow]++;
2497      }
2498    }
2499    CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
2500    int * indexRowL = coinFactorizationA_->indexRowL();
2501    int numberL = coinFactorizationA_->numberL();
2502    CoinBigIndex baseL = coinFactorizationA_->baseL();
2503    for (i=baseL;i<baseL+numberL;i++) {
2504      CoinBigIndex j;
2505      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
2506        int iRow = indexRowL[j];
2507        temp[iRow]++;
2508      }
2509    }
2510    for (i=0;i<numberRows;i++) {
2511      int number = temp[i];
2512      int iPermute = permuteBack[i];
2513      weights[iPermute]=number;
2514    }
2515    delete [] temp;
2516  } else {
2517    int i;
2518    for (i=0;i<numberRows;i++) {
2519      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
2520      //number = startRowL[i+1]-startRowL[i]+1;
2521      //number = numberInRow[i]+1;
2522      int iPermute = permuteBack[i];
2523      weights[iPermute]=number;
2524    }
2525  }
2526#ifdef CLP_FACTORIZATION_INSTRUMENT
2527  factorization_instrument(8);
2528#endif
2529}
2530#endif
Note: See TracBrowser for help on using the repository browser.