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

Last change on this file since 1958 was 1903, checked in by stefan, 7 years ago

code from John to check that fillBasis works correctly (need to be enabled via define)

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