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

Last change on this file since 2470 was 2385, checked in by unxusr, 9 months ago

formatting

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