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

Last change on this file since 1616 was 1616, checked in by forrest, 9 years ago

for partial skip of etas

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