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

Last change on this file since 2078 was 2078, checked in by forrest, 6 years ago

make faster on difficult problems

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