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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 121.0 KB
Line 
1/* $Id: ClpFactorization.cpp 1665 2011-01-04 17:55:54Z lou $ */
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     }
1402}
1403int
1404ClpFactorization::factorize ( ClpSimplex * model,
1405                              int solveType, bool valuesPass)
1406{
1407#ifdef CLP_REUSE_ETAS
1408     model_= model;
1409#endif
1410     //if ((model->specialOptions()&16384))
1411     //printf("factor at %d iterations\n",model->numberIterations());
1412     ClpMatrixBase * matrix = model->clpMatrix();
1413     int numberRows = model->numberRows();
1414     int numberColumns = model->numberColumns();
1415     if (!numberRows)
1416          return 0;
1417#ifdef CLP_FACTORIZATION_INSTRUMENT
1418     factorization_instrument(-1);
1419#endif
1420     bool anyChanged = false;
1421     if (coinFactorizationB_) {
1422          coinFactorizationB_->setStatus(-99);
1423          int * pivotVariable = model->pivotVariable();
1424          //returns 0 -okay, -1 singular, -2 too many in basis */
1425          // allow dense
1426          int solveMode = 2;
1427          if (model->numberIterations())
1428               solveMode += 8;
1429          if (valuesPass)
1430               solveMode += 4;
1431          coinFactorizationB_->setSolveMode(solveMode);
1432          while (status() < -98) {
1433
1434               int i;
1435               int numberBasic = 0;
1436               int numberRowBasic;
1437               // Move pivot variables across if they look good
1438               int * pivotTemp = model->rowArray(0)->getIndices();
1439               assert (!model->rowArray(0)->getNumElements());
1440               if (!matrix->rhsOffset(model)) {
1441                    // Seems to prefer things in order so quickest
1442                    // way is to go though like this
1443                    for (i = 0; i < numberRows; i++) {
1444                         if (model->getRowStatus(i) == ClpSimplex::basic)
1445                              pivotTemp[numberBasic++] = i;
1446                    }
1447                    numberRowBasic = numberBasic;
1448                    /* Put column basic variables into pivotVariable
1449                       This is done by ClpMatrixBase to allow override for gub
1450                    */
1451                    matrix->generalExpanded(model, 0, numberBasic);
1452               } else {
1453                    // Long matrix - do a different way
1454                    bool fullSearch = false;
1455                    for (i = 0; i < numberRows; i++) {
1456                         int iPivot = pivotVariable[i];
1457                         if (iPivot >= numberColumns) {
1458                              pivotTemp[numberBasic++] = iPivot - numberColumns;
1459                         }
1460                    }
1461                    numberRowBasic = numberBasic;
1462                    for (i = 0; i < numberRows; i++) {
1463                         int iPivot = pivotVariable[i];
1464                         if (iPivot < numberColumns) {
1465                              if (iPivot >= 0) {
1466                                   pivotTemp[numberBasic++] = iPivot;
1467                              } else {
1468                                   // not full basis
1469                                   fullSearch = true;
1470                                   break;
1471                              }
1472                         }
1473                    }
1474                    if (fullSearch) {
1475                         // do slow way
1476                         numberBasic = 0;
1477                         for (i = 0; i < numberRows; i++) {
1478                              if (model->getRowStatus(i) == ClpSimplex::basic)
1479                                   pivotTemp[numberBasic++] = i;
1480                         }
1481                         numberRowBasic = numberBasic;
1482                         /* Put column basic variables into pivotVariable
1483                            This is done by ClpMatrixBase to allow override for gub
1484                         */
1485                         matrix->generalExpanded(model, 0, numberBasic);
1486                    }
1487               }
1488               if (numberBasic > model->maximumBasic()) {
1489                    // Take out some
1490                    numberBasic = numberRowBasic;
1491                    for (int i = 0; i < numberColumns; i++) {
1492                         if (model->getColumnStatus(i) == ClpSimplex::basic) {
1493                              if (numberBasic < numberRows)
1494                                   numberBasic++;
1495                              else
1496                                   model->setColumnStatus(i, ClpSimplex::superBasic);
1497                         }
1498                    }
1499                    numberBasic = numberRowBasic;
1500                    matrix->generalExpanded(model, 0, numberBasic);
1501               } else if (numberBasic < numberRows) {
1502                    // add in some
1503                    int needed = numberRows - numberBasic;
1504                    // move up columns
1505                    for (i = numberBasic - 1; i >= numberRowBasic; i--)
1506                         pivotTemp[i+needed] = pivotTemp[i];
1507                    numberRowBasic = 0;
1508                    numberBasic = numberRows;
1509                    for (i = 0; i < numberRows; i++) {
1510                         if (model->getRowStatus(i) == ClpSimplex::basic) {
1511                              pivotTemp[numberRowBasic++] = i;
1512                         } else if (needed) {
1513                              needed--;
1514                              model->setRowStatus(i, ClpSimplex::basic);
1515                              pivotTemp[numberRowBasic++] = i;
1516                         }
1517                    }
1518               }
1519               CoinBigIndex numberElements = numberRowBasic;
1520
1521               // compute how much in basis
1522               // can change for gub
1523               int numberColumnBasic = numberBasic - numberRowBasic;
1524
1525               numberElements += matrix->countBasis(pivotTemp + numberRowBasic,
1526                                                    numberColumnBasic);
1527               // Not needed for dense
1528               numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1529               int numberIterations = model->numberIterations();
1530               coinFactorizationB_->setUsefulInformation(&numberIterations, 0);
1531               coinFactorizationB_->getAreas ( numberRows,
1532                                               numberRowBasic + numberColumnBasic, numberElements,
1533                                               2 * numberElements );
1534               // Fill in counts so we can skip part of preProcess
1535               // This is NOT needed for dense but would be needed for later versions
1536               CoinFactorizationDouble * elementU;
1537               int * indexRowU;
1538               CoinBigIndex * startColumnU;
1539               int * numberInRow;
1540               int * numberInColumn;
1541               elementU = coinFactorizationB_->elements();
1542               indexRowU = coinFactorizationB_->indices();
1543               startColumnU = coinFactorizationB_->starts();
1544#ifndef COIN_FAST_CODE
1545               double slackValue;
1546               slackValue = coinFactorizationB_->slackValue();
1547#else
1548#define slackValue -1.0
1549#endif
1550               numberInRow = coinFactorizationB_->numberInRow();
1551               numberInColumn = coinFactorizationB_->numberInColumn();
1552               CoinZeroN ( numberInRow, numberRows  );
1553               CoinZeroN ( numberInColumn, numberRows );
1554               for (i = 0; i < numberRowBasic; i++) {
1555                    int iRow = pivotTemp[i];
1556                    // Change pivotTemp to correct sequence
1557                    pivotTemp[i] = iRow + numberColumns;
1558                    indexRowU[i] = iRow;
1559                    startColumnU[i] = i;
1560                    elementU[i] = slackValue;
1561                    numberInRow[iRow] = 1;
1562                    numberInColumn[i] = 1;
1563               }
1564               startColumnU[numberRowBasic] = numberRowBasic;
1565               // can change for gub so redo
1566               numberColumnBasic = numberBasic - numberRowBasic;
1567               matrix->fillBasis(model,
1568                                 pivotTemp + numberRowBasic,
1569                                 numberColumnBasic,
1570                                 indexRowU,
1571                                 startColumnU + numberRowBasic,
1572                                 numberInRow,
1573                                 numberInColumn + numberRowBasic,
1574                                 elementU);
1575               // recompute number basic
1576               numberBasic = numberRowBasic + numberColumnBasic;
1577               for (i = numberBasic; i < numberRows; i++)
1578                    pivotTemp[i] = -1; // mark not there
1579               if (numberBasic)
1580                    numberElements = startColumnU[numberBasic-1]
1581                                     + numberInColumn[numberBasic-1];
1582               else
1583                    numberElements = 0;
1584               coinFactorizationB_->preProcess ( );
1585               coinFactorizationB_->factor (  );
1586               if (coinFactorizationB_->status() == -1 &&
1587                         (coinFactorizationB_->solveMode() % 3) != 0) {
1588                    int solveMode = coinFactorizationB_->solveMode();
1589                    solveMode -= solveMode % 3; // so bottom will be 0
1590                    coinFactorizationB_->setSolveMode(solveMode);
1591                    coinFactorizationB_->setStatus(-99);
1592               }
1593               if (coinFactorizationB_->status() == -99)
1594                    continue;
1595               // If we get here status is 0 or -1
1596               if (coinFactorizationB_->status() == 0 && numberBasic == numberRows) {
1597                    coinFactorizationB_->postProcess(pivotTemp, pivotVariable);
1598               } else if (solveType == 0 || solveType == 2/*||solveType==1*/) {
1599                    // Change pivotTemp to be correct list
1600                    anyChanged = true;
1601                    coinFactorizationB_->makeNonSingular(pivotTemp, numberColumns);
1602                    double * columnLower = model->lowerRegion();
1603                    double * columnUpper = model->upperRegion();
1604                    double * columnActivity = model->solutionRegion();
1605                    double * rowLower = model->lowerRegion(0);
1606                    double * rowUpper = model->upperRegion(0);
1607                    double * rowActivity = model->solutionRegion(0);
1608                    //redo basis - first take ALL out
1609                    int iColumn;
1610                    double largeValue = model->largeValue();
1611                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1612                         if (model->getColumnStatus(iColumn) == ClpSimplex::basic) {
1613                              // take out
1614                              if (!valuesPass) {
1615                                   double lower = columnLower[iColumn];
1616                                   double upper = columnUpper[iColumn];
1617                                   double value = columnActivity[iColumn];
1618                                   if (lower > -largeValue || upper < largeValue) {
1619                                        if (fabs(value - lower) < fabs(value - upper)) {
1620                                             model->setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1621                                             columnActivity[iColumn] = lower;
1622                                        } else {
1623                                             model->setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1624                                             columnActivity[iColumn] = upper;
1625                                        }
1626                                   } else {
1627                                        model->setColumnStatus(iColumn, ClpSimplex::isFree);
1628                                   }
1629                              } else {
1630                                   model->setColumnStatus(iColumn, ClpSimplex::superBasic);
1631                              }
1632                         }
1633                    }
1634                    int iRow;
1635                    for (iRow = 0; iRow < numberRows; iRow++) {
1636                         if (model->getRowStatus(iRow) == ClpSimplex::basic) {
1637                              // take out
1638                              if (!valuesPass) {
1639                                   double lower = columnLower[iRow];
1640                                   double upper = columnUpper[iRow];
1641                                   double value = columnActivity[iRow];
1642                                   if (lower > -largeValue || upper < largeValue) {
1643                                        if (fabs(value - lower) < fabs(value - upper)) {
1644                                             model->setRowStatus(iRow, ClpSimplex::atLowerBound);
1645                                             columnActivity[iRow] = lower;
1646                                        } else {
1647                                             model->setRowStatus(iRow, ClpSimplex::atUpperBound);
1648                                             columnActivity[iRow] = upper;
1649                                        }
1650                                   } else {
1651                                        model->setRowStatus(iRow, ClpSimplex::isFree);
1652                                   }
1653                              } else {
1654                                   model->setRowStatus(iRow, ClpSimplex::superBasic);
1655                              }
1656                         }
1657                    }
1658                    for (iRow = 0; iRow < numberRows; iRow++) {
1659                         int iSequence = pivotTemp[iRow];
1660                         assert (iSequence >= 0);
1661                         // basic
1662                         model->setColumnStatus(iSequence, ClpSimplex::basic);
1663                    }
1664                    // signal repeat
1665                    coinFactorizationB_->setStatus(-99);
1666                    // set fixed if they are
1667                    for (iRow = 0; iRow < numberRows; iRow++) {
1668                         if (model->getRowStatus(iRow) != ClpSimplex::basic ) {
1669                              if (rowLower[iRow] == rowUpper[iRow]) {
1670                                   rowActivity[iRow] = rowLower[iRow];
1671                                   model->setRowStatus(iRow, ClpSimplex::isFixed);
1672                              }
1673                         }
1674                    }
1675                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1676                         if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) {
1677                              if (columnLower[iColumn] == columnUpper[iColumn]) {
1678                                   columnActivity[iColumn] = columnLower[iColumn];
1679                                   model->setColumnStatus(iColumn, ClpSimplex::isFixed);
1680                              }
1681                         }
1682                    }
1683               }
1684          }
1685#ifdef CLP_DEBUG
1686          // check basic
1687          CoinIndexedVector region1(2 * numberRows);
1688          CoinIndexedVector region2B(2 * numberRows);
1689          int iPivot;
1690          double * arrayB = region2B.denseVector();
1691          int i;
1692          for (iPivot = 0; iPivot < numberRows; iPivot++) {
1693               int iSequence = pivotVariable[iPivot];
1694               model->unpack(&region2B, iSequence);
1695               coinFactorizationB_->updateColumn(&region1, &region2B);
1696               if (fabs(arrayB[iPivot] - 1.0) < 1.0e-4) {
1697                    // OK?
1698                    arrayB[iPivot] = 0.0;
1699               } else {
1700                    assert (fabs(arrayB[iPivot]) < 1.0e-4);
1701                    for (i = 0; i < numberRows; i++) {
1702                         if (fabs(arrayB[i] - 1.0) < 1.0e-4)
1703                              break;
1704                    }
1705                    assert (i < numberRows);
1706                    printf("variable on row %d landed up on row %d\n", iPivot, i);
1707                    arrayB[i] = 0.0;
1708               }
1709               for (i = 0; i < numberRows; i++)
1710                    assert (fabs(arrayB[i]) < 1.0e-4);
1711               region2B.clear();
1712          }
1713#endif
1714#ifdef CLP_FACTORIZATION_INSTRUMENT
1715          factorization_instrument(2);
1716#endif
1717          if ( anyChanged && model->algorithm() < 0 && solveType > 0) {
1718               double dummyCost;
1719               static_cast<ClpSimplexDual *> (model)->changeBounds(3,
1720                         NULL, dummyCost);
1721          }
1722          return coinFactorizationB_->status();
1723     }
1724     // If too many compressions increase area
1725     if (coinFactorizationA_->pivots() > 1 && coinFactorizationA_->numberCompressions() * 10 > coinFactorizationA_->pivots() + 10) {
1726          coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor() * 1.1);
1727     }
1728     //int numberPivots=coinFactorizationA_->pivots();
1729#if 0
1730     if (model->algorithm() > 0)
1731          numberSave = -1;
1732#endif
1733#ifndef SLIM_CLP
1734     if (!networkBasis_ || doCheck) {
1735#endif
1736          coinFactorizationA_->setStatus(-99);
1737          int * pivotVariable = model->pivotVariable();
1738          int nTimesRound = 0;
1739          //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
1740          while (coinFactorizationA_->status() < -98) {
1741               nTimesRound++;
1742
1743               int i;
1744               int numberBasic = 0;
1745               int numberRowBasic;
1746               // Move pivot variables across if they look good
1747               int * pivotTemp = model->rowArray(0)->getIndices();
1748               assert (!model->rowArray(0)->getNumElements());
1749               if (!matrix->rhsOffset(model)) {
1750#if 0
1751                    if (numberSave > 0) {
1752                         int nStill = 0;
1753                         int nAtBound = 0;
1754                         int nZeroDual = 0;
1755                         CoinIndexedVector * array = model->rowArray(3);
1756                         CoinIndexedVector * objArray = model->columnArray(1);
1757                         array->clear();
1758                         objArray->clear();
1759                         double * cost = model->costRegion();
1760                         double tolerance = model->primalTolerance();
1761                         double offset = 0.0;
1762                         for (i = 0; i < numberRows; i++) {
1763                              int iPivot = pivotVariable[i];
1764                              if (iPivot < numberColumns && isDense(iPivot)) {
1765                                   if (model->getColumnStatus(iPivot) == ClpSimplex::basic) {
1766                                        nStill++;
1767                                        double value = model->solutionRegion()[iPivot];
1768                                        double dual = model->dualRowSolution()[i];
1769                                        double lower = model->lowerRegion()[iPivot];
1770                                        double upper = model->upperRegion()[iPivot];
1771                                        ClpSimplex::Status status;
1772                                        if (fabs(value - lower) < tolerance) {
1773                                             status = ClpSimplex::atLowerBound;
1774                                             nAtBound++;
1775                                        } else if (fabs(value - upper) < tolerance) {
1776                                             nAtBound++;
1777                                             status = ClpSimplex::atUpperBound;
1778                                        } else if (value > lower && value < upper) {
1779                                             status = ClpSimplex::superBasic;
1780                                        } else {
1781                                             status = ClpSimplex::basic;
1782                                        }
1783                                        if (status != ClpSimplex::basic) {
1784                                             if (model->getRowStatus(i) != ClpSimplex::basic) {
1785                                                  model->setColumnStatus(iPivot, ClpSimplex::atLowerBound);
1786                                                  model->setRowStatus(i, ClpSimplex::basic);
1787                                                  pivotVariable[i] = i + numberColumns;
1788                                                  model->dualRowSolution()[i] = 0.0;
1789                                                  model->djRegion(0)[i] = 0.0;
1790                                                  array->add(i, dual);
1791                                                  offset += dual * model->solutionRegion(0)[i];
1792                                             }
1793                                        }
1794                                        if (fabs(dual) < 1.0e-5)
1795                                             nZeroDual++;
1796                                   }
1797                              }
1798                         }
1799                         printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
1800                                numberSave, nStill, nAtBound, nZeroDual, offset);
1801                         if (array->getNumElements()) {
1802                              // modify costs
1803                              model->clpMatrix()->transposeTimes(model, 1.0, array, model->columnArray(0),
1804                                                                 objArray);
1805                              array->clear();
1806                              int n = objArray->getNumElements();
1807                              int * indices = objArray->getIndices();
1808                              double * elements = objArray->denseVector();
1809                              for (i = 0; i < n; i++) {
1810                                   int iColumn = indices[i];
1811                                   cost[iColumn] -= elements[iColumn];
1812                                   elements[iColumn] = 0.0;
1813                              }
1814                              objArray->setNumElements(0);
1815                         }
1816                    }
1817#endif
1818                    // Seems to prefer things in order so quickest
1819                    // way is to go though like this
1820                    for (i = 0; i < numberRows; i++) {
1821                         if (model->getRowStatus(i) == ClpSimplex::basic)
1822                              pivotTemp[numberBasic++] = i;
1823                    }
1824                    numberRowBasic = numberBasic;
1825                    /* Put column basic variables into pivotVariable
1826                       This is done by ClpMatrixBase to allow override for gub
1827                    */
1828                    matrix->generalExpanded(model, 0, numberBasic);
1829               } else {
1830                    // Long matrix - do a different way
1831                    bool fullSearch = false;
1832                    for (i = 0; i < numberRows; i++) {
1833                         int iPivot = pivotVariable[i];
1834                         if (iPivot >= numberColumns) {
1835                              pivotTemp[numberBasic++] = iPivot - numberColumns;
1836                         }
1837                    }
1838                    numberRowBasic = numberBasic;
1839                    for (i = 0; i < numberRows; i++) {
1840                         int iPivot = pivotVariable[i];
1841                         if (iPivot < numberColumns) {
1842                              if (iPivot >= 0) {
1843                                   pivotTemp[numberBasic++] = iPivot;
1844                              } else {
1845                                   // not full basis
1846                                   fullSearch = true;
1847                                   break;
1848                              }
1849                         }
1850                    }
1851                    if (fullSearch) {
1852                         // do slow way
1853                         numberBasic = 0;
1854                         for (i = 0; i < numberRows; i++) {
1855                              if (model->getRowStatus(i) == ClpSimplex::basic)
1856                                   pivotTemp[numberBasic++] = i;
1857                         }
1858                         numberRowBasic = numberBasic;
1859                         /* Put column basic variables into pivotVariable
1860                            This is done by ClpMatrixBase to allow override for gub
1861                         */
1862                         matrix->generalExpanded(model, 0, numberBasic);
1863                    }
1864               }
1865               if (numberBasic > model->maximumBasic()) {
1866#if 0 // ndef NDEBUG
1867                    printf("%d basic - should only be %d\n",
1868                           numberBasic, numberRows);
1869#endif
1870                    // Take out some
1871                    numberBasic = numberRowBasic;
1872                    for (int i = 0; i < numberColumns; i++) {
1873                         if (model->getColumnStatus(i) == ClpSimplex::basic) {
1874                              if (numberBasic < numberRows)
1875                                   numberBasic++;
1876                              else
1877                                   model->setColumnStatus(i, ClpSimplex::superBasic);
1878                         }
1879                    }
1880                    numberBasic = numberRowBasic;
1881                    matrix->generalExpanded(model, 0, numberBasic);
1882               }
1883#ifndef SLIM_CLP
1884               // see if matrix a network
1885#ifndef NO_RTTI
1886               ClpNetworkMatrix* networkMatrix =
1887                    dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
1888#else
1889ClpNetworkMatrix* networkMatrix = NULL;
1890if (model->clpMatrix()->type() == 11)
1891     networkMatrix =
1892          static_cast< ClpNetworkMatrix*>(model->clpMatrix());
1893#endif
1894               // If network - still allow ordinary factorization first time for laziness
1895               if (networkMatrix)
1896                    coinFactorizationA_->setBiasLU(0); // All to U if network
1897               //int saveMaximumPivots = maximumPivots();
1898               delete networkBasis_;
1899               networkBasis_ = NULL;
1900               if (networkMatrix && !doCheck)
1901                    maximumPivots(1);
1902#endif
1903               //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
1904               while (coinFactorizationA_->status() == -99) {
1905                    // maybe for speed will be better to leave as many regions as possible
1906                    coinFactorizationA_->gutsOfDestructor();
1907                    coinFactorizationA_->gutsOfInitialize(2);
1908                    CoinBigIndex numberElements = numberRowBasic;
1909
1910                    // compute how much in basis
1911
1912                    int i;
1913                    // can change for gub
1914                    int numberColumnBasic = numberBasic - numberRowBasic;
1915
1916                    numberElements += matrix->countBasis( pivotTemp + numberRowBasic,
1917                                                          numberColumnBasic);
1918                    // and recompute as network side say different
1919                    if (model->numberIterations())
1920                         numberRowBasic = numberBasic - numberColumnBasic;
1921                    numberElements = 3 * numberBasic + 3 * numberElements + 20000;
1922                    coinFactorizationA_->getAreas ( numberRows,
1923                                                    numberRowBasic + numberColumnBasic, numberElements,
1924                                                    2 * numberElements );
1925                    //fill
1926                    // Fill in counts so we can skip part of preProcess
1927                    int * numberInRow = coinFactorizationA_->numberInRow();
1928                    int * numberInColumn = coinFactorizationA_->numberInColumn();
1929                    CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
1930                    CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
1931                    CoinFactorizationDouble * elementU = coinFactorizationA_->elementU();
1932                    int * indexRowU = coinFactorizationA_->indexRowU();
1933                    CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
1934#ifndef COIN_FAST_CODE
1935                    double slackValue = coinFactorizationA_->slackValue();
1936#endif
1937                    for (i = 0; i < numberRowBasic; i++) {
1938                         int iRow = pivotTemp[i];
1939                         indexRowU[i] = iRow;
1940                         startColumnU[i] = i;
1941                         elementU[i] = slackValue;
1942                         numberInRow[iRow] = 1;
1943                         numberInColumn[i] = 1;
1944                    }
1945                    startColumnU[numberRowBasic] = numberRowBasic;
1946                    // can change for gub so redo
1947                    numberColumnBasic = numberBasic - numberRowBasic;
1948                    matrix->fillBasis(model,
1949                                      pivotTemp + numberRowBasic,
1950                                      numberColumnBasic,
1951                                      indexRowU,
1952                                      startColumnU + numberRowBasic,
1953                                      numberInRow,
1954                                      numberInColumn + numberRowBasic,
1955                                      elementU);
1956#if 0
1957                    {
1958                         printf("%d row basic, %d column basic\n", numberRowBasic, numberColumnBasic);
1959                         for (int i = 0; i < numberElements; i++)
1960                              printf("row %d col %d value %g\n", indexRowU[i], indexColumnU_[i],
1961                                     elementU[i]);
1962                    }
1963#endif
1964                    // recompute number basic
1965                    numberBasic = numberRowBasic + numberColumnBasic;
1966                    if (numberBasic)
1967                         numberElements = startColumnU[numberBasic-1]
1968                                          + numberInColumn[numberBasic-1];
1969                    else
1970                         numberElements = 0;
1971                    coinFactorizationA_->setNumberElementsU(numberElements);
1972                    //saveFactorization("dump.d");
1973                    if (coinFactorizationA_->biasLU() >= 3 || coinFactorizationA_->numberRows() != coinFactorizationA_->numberColumns())
1974                         coinFactorizationA_->preProcess ( 2 );
1975                    else
1976                         coinFactorizationA_->preProcess ( 3 ); // no row copy
1977                    coinFactorizationA_->factor (  );
1978                    if (coinFactorizationA_->status() == -99) {
1979                         // get more memory
1980                         coinFactorizationA_->areaFactor(2.0 * coinFactorizationA_->areaFactor());
1981                    } else if (coinFactorizationA_->status() == -1 &&
1982                               (model->numberIterations() == 0 || nTimesRound > 2) &&
1983                               coinFactorizationA_->denseThreshold()) {
1984                         // Round again without dense
1985                         coinFactorizationA_->setDenseThreshold(0);
1986                         coinFactorizationA_->setStatus(-99);
1987                    }
1988               }
1989               // If we get here status is 0 or -1
1990               if (coinFactorizationA_->status() == 0) {
1991                    // We may need to tamper with order and redo - e.g. network with side
1992                    int useNumberRows = numberRows;
1993                    // **** we will also need to add test in dual steepest to do
1994                    // as we do for network
1995                    matrix->generalExpanded(model, 12, useNumberRows);
1996                    const int * permuteBack = coinFactorizationA_->permuteBack();
1997                    const int * back = coinFactorizationA_->pivotColumnBack();
1998                    //int * pivotTemp = pivotColumn_.array();
1999                    //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
2000#ifndef NDEBUG
2001                    CoinFillN(pivotVariable, numberRows, -1);
2002#endif
2003                    // Redo pivot order
2004                    for (i = 0; i < numberRowBasic; i++) {
2005                         int k = pivotTemp[i];
2006                         // so rowIsBasic[k] would be permuteBack[back[i]]
2007                         int j = permuteBack[back[i]];
2008                         assert (pivotVariable[j] == -1);
2009                         pivotVariable[j] = k + numberColumns;
2010                    }
2011                    for (; i < useNumberRows; i++) {
2012                         int k = pivotTemp[i];
2013                         // so rowIsBasic[k] would be permuteBack[back[i]]
2014                         int j = permuteBack[back[i]];
2015                         assert (pivotVariable[j] == -1);
2016                         pivotVariable[j] = k;
2017                    }
2018#if 0
2019                    if (numberSave >= 0) {
2020                         numberSave = numberDense_;
2021                         memset(saveList, 0, ((coinFactorizationA_->numberRows() + 31) >> 5)*sizeof(int));
2022                         for (i = coinFactorizationA_->numberRows() - numberSave; i < coinFactorizationA_->numberRows(); i++) {
2023                              int k = pivotTemp[pivotColumn_.array()[i]];
2024                              setDense(k);
2025                         }
2026                    }
2027#endif
2028                    // Set up permutation vector
2029                    // these arrays start off as copies of permute
2030                    // (and we could use permute_ instead of pivotColumn (not back though))
2031                    ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn()  );
2032                    ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack()  );
2033#ifndef SLIM_CLP
2034                    if (networkMatrix) {
2035                         maximumPivots(CoinMax(2000, maximumPivots()));
2036                         // redo arrays
2037                         for (int iRow = 0; iRow < 4; iRow++) {
2038                              int length = model->numberRows() + maximumPivots();
2039                              if (iRow == 3 || model->objectiveAsObject()->type() > 1)
2040                                   length += model->numberColumns();
2041                              model->rowArray(iRow)->reserve(length);
2042                         }
2043                         // create network factorization
2044                         if (doCheck)
2045                              delete networkBasis_; // temp
2046                         networkBasis_ = new ClpNetworkBasis(model, coinFactorizationA_->numberRows(),
2047                                                             coinFactorizationA_->pivotRegion(),
2048                                                             coinFactorizationA_->permuteBack(),
2049                                                             coinFactorizationA_->startColumnU(),
2050                                                             coinFactorizationA_->numberInColumn(),
2051                                                             coinFactorizationA_->indexRowU(),
2052                                                             coinFactorizationA_->elementU());
2053                         // kill off arrays in ordinary factorization
2054                         if (!doCheck) {
2055                              coinFactorizationA_->gutsOfDestructor();
2056                              // but make sure coinFactorizationA_->numberRows() set
2057                              coinFactorizationA_->setNumberRows(model->numberRows());
2058                              coinFactorizationA_->setStatus(0);
2059#if 0
2060                              // but put back permute arrays so odd things will work
2061                              int numberRows = model->numberRows();
2062                              pivotColumnBack_ = new int [numberRows];
2063                              permute_ = new int [numberRows];
2064                              int i;
2065                              for (i = 0; i < numberRows; i++) {
2066                                   pivotColumnBack_[i] = i;
2067                                   permute_[i] = i;
2068                              }
2069#endif
2070                         }
2071                    } else {
2072#endif
2073                         // See if worth going sparse and when
2074                         coinFactorizationA_->checkSparse();
2075#ifndef SLIM_CLP
2076                    }
2077#endif
2078               } else if (coinFactorizationA_->status() == -1 && (solveType == 0 || solveType == 2)) {
2079                    // This needs redoing as it was merged coding - does not need array
2080                    int numberTotal = numberRows + numberColumns;
2081                    int * isBasic = new int [numberTotal];
2082                    int * rowIsBasic = isBasic + numberColumns;
2083                    int * columnIsBasic = isBasic;
2084                    for (i = 0; i < numberTotal; i++)
2085                         isBasic[i] = -1;
2086                    for (i = 0; i < numberRowBasic; i++) {
2087                         int iRow = pivotTemp[i];
2088                         rowIsBasic[iRow] = 1;
2089                    }
2090                    for (; i < numberBasic; i++) {
2091                         int iColumn = pivotTemp[i];
2092                         columnIsBasic[iColumn] = 1;
2093                    }
2094                    numberBasic = 0;
2095                    for (i = 0; i < numberRows; i++)
2096                         pivotVariable[i] = -1;
2097                    // mark as basic or non basic
2098                    const int * pivotColumn = coinFactorizationA_->pivotColumn();
2099                    for (i = 0; i < numberRows; i++) {
2100                         if (rowIsBasic[i] >= 0) {
2101                              if (pivotColumn[numberBasic] >= 0) {
2102                                   rowIsBasic[i] = pivotColumn[numberBasic];
2103                              } else {
2104                                   rowIsBasic[i] = -1;
2105                                   model->setRowStatus(i, ClpSimplex::superBasic);
2106                              }
2107                              numberBasic++;
2108                         }
2109                    }
2110                    for (i = 0; i < numberColumns; i++) {
2111                         if (columnIsBasic[i] >= 0) {
2112                              if (pivotColumn[numberBasic] >= 0)
2113                                   columnIsBasic[i] = pivotColumn[numberBasic];
2114                              else
2115                                   columnIsBasic[i] = -1;
2116                              numberBasic++;
2117                         }
2118                    }
2119                    // leave pivotVariable in useful form for cleaning basis
2120                    int * pivotVariable = model->pivotVariable();
2121                    for (i = 0; i < numberRows; i++) {
2122                         pivotVariable[i] = -1;
2123                    }
2124
2125                    for (i = 0; i < numberRows; i++) {
2126                         if (model->getRowStatus(i) == ClpSimplex::basic) {
2127                              int iPivot = rowIsBasic[i];
2128                              if (iPivot >= 0)
2129                                   pivotVariable[iPivot] = i + numberColumns;
2130                         }
2131                    }
2132                    for (i = 0; i < numberColumns; i++) {
2133                         if (model->getColumnStatus(i) == ClpSimplex::basic) {
2134                              int iPivot = columnIsBasic[i];
2135                              if (iPivot >= 0)
2136                                   pivotVariable[iPivot] = i;
2137                         }
2138                    }
2139                    delete [] isBasic;
2140                    double * columnLower = model->lowerRegion();
2141                    double * columnUpper = model->upperRegion();
2142                    double * columnActivity = model->solutionRegion();
2143                    double * rowLower = model->lowerRegion(0);
2144                    double * rowUpper = model->upperRegion(0);
2145                    double * rowActivity = model->solutionRegion(0);
2146                    //redo basis - first take ALL columns out
2147                    int iColumn;
2148                    double largeValue = model->largeValue();
2149                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2150                         if (model->getColumnStatus(iColumn) == ClpSimplex::basic) {
2151                              // take out
2152                              if (!valuesPass) {
2153                                   double lower = columnLower[iColumn];
2154                                   double upper = columnUpper[iColumn];
2155                                   double value = columnActivity[iColumn];
2156                                   if (lower > -largeValue || upper < largeValue) {
2157                                        if (fabs(value - lower) < fabs(value - upper)) {
2158                                             model->setColumnStatus(iColumn, ClpSimplex::atLowerBound);
2159                                             columnActivity[iColumn] = lower;
2160                                        } else {
2161                                             model->setColumnStatus(iColumn, ClpSimplex::atUpperBound);
2162                                             columnActivity[iColumn] = upper;
2163                                        }
2164                                   } else {
2165                                        model->setColumnStatus(iColumn, ClpSimplex::isFree);
2166                                   }
2167                              } else {
2168                                   model->setColumnStatus(iColumn, ClpSimplex::superBasic);
2169                              }
2170                         }
2171                    }
2172                    int iRow;
2173                    for (iRow = 0; iRow < numberRows; iRow++) {
2174                         int iSequence = pivotVariable[iRow];
2175                         if (iSequence >= 0) {
2176                              // basic
2177                              if (iSequence >= numberColumns) {
2178                                   // slack in - leave
2179                                   //assert (iSequence-numberColumns==iRow);
2180                              } else {
2181                                   assert(model->getRowStatus(iRow) != ClpSimplex::basic);
2182                                   // put back structural
2183                                   model->setColumnStatus(iSequence, ClpSimplex::basic);
2184                              }
2185                         } else {
2186                              // put in slack
2187                              model->setRowStatus(iRow, ClpSimplex::basic);
2188                         }
2189                    }
2190                    // Put back any key variables for gub
2191                    int dummy;
2192                    matrix->generalExpanded(model, 1, dummy);
2193                    // signal repeat
2194                    coinFactorizationA_->setStatus(-99);
2195                    // set fixed if they are
2196                    for (iRow = 0; iRow < numberRows; iRow++) {
2197                         if (model->getRowStatus(iRow) != ClpSimplex::basic ) {
2198                              if (rowLower[iRow] == rowUpper[iRow]) {
2199                                   rowActivity[iRow] = rowLower[iRow];
2200                                   model->setRowStatus(iRow, ClpSimplex::isFixed);
2201                              }
2202                         }
2203                    }
2204                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2205                         if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) {
2206                              if (columnLower[iColumn] == columnUpper[iColumn]) {
2207                                   columnActivity[iColumn] = columnLower[iColumn];
2208                                   model->setColumnStatus(iColumn, ClpSimplex::isFixed);
2209                              }
2210                         }
2211                    }
2212               }
2213          }
2214#ifndef SLIM_CLP
2215     } else {
2216          // network - fake factorization - do nothing
2217          coinFactorizationA_->setStatus(0);
2218          coinFactorizationA_->setPivots(0);
2219     }
2220#endif
2221#ifndef SLIM_CLP
2222     if (!coinFactorizationA_->status()) {
2223          // take out part if quadratic
2224          if (model->algorithm() == 2) {
2225               ClpObjective * obj = model->objectiveAsObject();
2226#ifndef NDEBUG
2227               ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
2228               assert (quadraticObj);
2229#else
2230ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
2231#endif
2232               CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2233               int numberXColumns = quadratic->getNumCols();
2234               assert (numberXColumns < numberColumns);
2235               int base = numberColumns - numberXColumns;
2236               int * which = new int [numberXColumns];
2237               int * pivotVariable = model->pivotVariable();
2238               int * permute = pivotColumn();
2239               int i;
2240               int n = 0;
2241               for (i = 0; i < numberRows; i++) {
2242                    int iSj = pivotVariable[i] - base;
2243                    if (iSj >= 0 && iSj < numberXColumns)
2244                         which[n++] = permute[i];
2245               }
2246               if (n)
2247                    coinFactorizationA_->emptyRows(n, which);
2248               delete [] which;
2249          }
2250     }
2251#endif
2252#ifdef CLP_FACTORIZATION_INSTRUMENT
2253     factorization_instrument(2);
2254#endif
2255     return coinFactorizationA_->status();
2256}
2257/* Replaces one Column in basis,
2258   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
2259   If checkBeforeModifying is true will do all accuracy checks
2260   before modifying factorization.  Whether to set this depends on
2261   speed considerations.  You could just do this on first iteration
2262   after factorization and thereafter re-factorize
2263   partial update already in U */
2264int
2265ClpFactorization::replaceColumn ( const ClpSimplex * model,
2266                                  CoinIndexedVector * regionSparse,
2267                                  CoinIndexedVector * tableauColumn,
2268                                  int pivotRow,
2269                                  double pivotCheck ,
2270                                  bool checkBeforeModifying,
2271                                  double acceptablePivot)
2272{
2273#ifndef SLIM_CLP
2274     if (!networkBasis_) {
2275#endif
2276#ifdef CLP_FACTORIZATION_INSTRUMENT
2277          factorization_instrument(-1);
2278#endif
2279          int returnCode;
2280          // see if FT
2281          if (!coinFactorizationA_ || coinFactorizationA_->forrestTomlin()) {
2282               if (coinFactorizationA_) {
2283                    returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2284                                 pivotRow,
2285                                 pivotCheck,
2286                                 checkBeforeModifying,
2287                                 acceptablePivot);
2288               } else {
2289                    bool tab = coinFactorizationB_->wantsTableauColumn();
2290#ifdef CLP_REUSE_ETAS
2291                    int tempInfo[2];
2292                    tempInfo[1] = model_->sequenceOut();
2293#else
2294                    int tempInfo[1];
2295#endif
2296                    tempInfo[0] = model->numberIterations();
2297                    coinFactorizationB_->setUsefulInformation(tempInfo, 1);
2298                    returnCode =
2299                         coinFactorizationB_->replaceColumn(tab ? tableauColumn : regionSparse,
2300                                                            pivotRow,
2301                                                            pivotCheck,
2302                                                            checkBeforeModifying,
2303                                                            acceptablePivot);
2304#ifdef CLP_DEBUG
2305                    // check basic
2306                    int numberRows = coinFactorizationB_->numberRows();
2307                    CoinIndexedVector region1(2 * numberRows);
2308                    CoinIndexedVector region2A(2 * numberRows);
2309                    CoinIndexedVector region2B(2 * numberRows);
2310                    int iPivot;
2311                    double * arrayB = region2B.denseVector();
2312                    int * pivotVariable = model->pivotVariable();
2313                    int i;
2314                    for (iPivot = 0; iPivot < numberRows; iPivot++) {
2315                         int iSequence = pivotVariable[iPivot];
2316                         if (iPivot == pivotRow)
2317                              iSequence = model->sequenceIn();
2318                         model->unpack(&region2B, iSequence);
2319                         coinFactorizationB_->updateColumn(&region1, &region2B);
2320                         assert (fabs(arrayB[iPivot] - 1.0) < 1.0e-4);
2321                         arrayB[iPivot] = 0.0;
2322                         for (i = 0; i < numberRows; i++)
2323                              assert (fabs(arrayB[i]) < 1.0e-4);
2324                         region2B.clear();
2325                    }
2326#endif
2327               }
2328          } else {
2329               returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn,
2330                            pivotRow, pivotCheck); // Note array
2331          }
2332#ifdef CLP_FACTORIZATION_INSTRUMENT
2333          factorization_instrument(3);
2334#endif
2335          return returnCode;
2336
2337#ifndef SLIM_CLP
2338     } else {
2339          if (doCheck) {
2340               int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
2341                                pivotRow,
2342                                pivotCheck,
2343                                checkBeforeModifying,
2344                                acceptablePivot);
2345               networkBasis_->replaceColumn(regionSparse,
2346                                            pivotRow);
2347               return returnCode;
2348          } else {
2349               // increase number of pivots
2350               coinFactorizationA_->setPivots(coinFactorizationA_->pivots() + 1);
2351               return networkBasis_->replaceColumn(regionSparse,
2352                                                   pivotRow);
2353          }
2354     }
2355#endif
2356}
2357
2358/* Updates one column (FTRAN) from region2
2359   number returned is negative if no room
2360   region1 starts as zero and is zero at end */
2361int
2362ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
2363                                   CoinIndexedVector * regionSparse2)
2364{
2365#ifdef CLP_DEBUG
2366     regionSparse->checkClear();
2367#endif
2368     if (!numberRows())
2369          return 0;
2370#ifndef SLIM_CLP
2371     if (!networkBasis_) {
2372#endif
2373#ifdef CLP_FACTORIZATION_INSTRUMENT
2374          factorization_instrument(-1);
2375#endif
2376          int returnCode;
2377          if (coinFactorizationA_) {
2378               coinFactorizationA_->setCollectStatistics(true);
2379               returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2380                            regionSparse2);
2381               coinFactorizationA_->setCollectStatistics(false);
2382          } else {
2383#ifdef CLP_REUSE_ETAS
2384              int tempInfo[2];
2385              tempInfo[0] = model_->numberIterations();
2386              tempInfo[1] = model_->sequenceIn();
2387              coinFactorizationB_->setUsefulInformation(tempInfo, 2);
2388#endif
2389              returnCode = coinFactorizationB_->updateColumnFT(regionSparse,
2390                            regionSparse2);
2391          }
2392#ifdef CLP_FACTORIZATION_INSTRUMENT
2393          factorization_instrument(4);
2394#endif
2395          return returnCode;
2396#ifndef SLIM_CLP
2397     } else {
2398#ifdef CHECK_NETWORK
2399          CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2400          double * check = new double[coinFactorizationA_->numberRows()];
2401          int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
2402                           regionSparse2);
2403          networkBasis_->updateColumn(regionSparse, save, -1);
2404          int i;
2405          double * array = regionSparse2->denseVector();
2406          int * indices = regionSparse2->getIndices();
2407          int n = regionSparse2->getNumElements();
2408          memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double));
2409          double * array2 = save->denseVector();
2410          int * indices2 = save->getIndices();
2411          int n2 = save->getNumElements();
2412          assert (n == n2);
2413          if (save->packedMode()) {
2414               for (i = 0; i < n; i++) {
2415                    check[indices[i]] = array[i];
2416               }
2417               for (i = 0; i < n; i++) {
2418                    double value2 = array2[i];
2419                    assert (check[indices2[i]] == value2);
2420               }
2421          } else {
2422               int numberRows = coinFactorizationA_->numberRows();
2423               for (i = 0; i < numberRows; i++) {
2424                    double value1 = array[i];
2425                    double value2 = array2[i];
2426                    assert (value1 == value2);
2427               }
2428          }
2429          delete save;
2430          delete [] check;
2431          return returnCode;
2432#else
2433networkBasis_->updateColumn(regionSparse, regionSparse2, -1);
2434return 1;
2435#endif
2436     }
2437#endif
2438}
2439/* Updates one column (FTRAN) from region2
2440   number returned is negative if no room
2441   region1 starts as zero and is zero at end */
2442int
2443ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
2444                                 CoinIndexedVector * regionSparse2,
2445                                 bool noPermute) const
2446{
2447#ifdef CLP_DEBUG
2448     if (!noPermute)
2449          regionSparse->checkClear();
2450#endif
2451     if (!numberRows())
2452          return 0;
2453#ifndef SLIM_CLP
2454     if (!networkBasis_) {
2455#endif
2456#ifdef CLP_FACTORIZATION_INSTRUMENT
2457          factorization_instrument(-1);
2458#endif
2459          int returnCode;
2460          if (coinFactorizationA_) {
2461               coinFactorizationA_->setCollectStatistics(true);
2462               returnCode = coinFactorizationA_->updateColumn(regionSparse,
2463                            regionSparse2,
2464                            noPermute);
2465               coinFactorizationA_->setCollectStatistics(false);
2466          } else {
2467               returnCode = coinFactorizationB_->updateColumn(regionSparse,
2468                            regionSparse2,
2469                            noPermute);
2470          }
2471#ifdef CLP_FACTORIZATION_INSTRUMENT
2472          factorization_instrument(5);
2473#endif
2474          //#define PRINT_VECTOR
2475#ifdef PRINT_VECTOR
2476          printf("Update\n");
2477          regionSparse2->print();
2478#endif
2479          return returnCode;
2480#ifndef SLIM_CLP
2481     } else {
2482#ifdef CHECK_NETWORK
2483          CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2484          double * check = new double[coinFactorizationA_->numberRows()];
2485          int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2486                           regionSparse2,
2487                           noPermute);
2488          networkBasis_->updateColumn(regionSparse, save, -1);
2489          int i;
2490          double * array = regionSparse2->denseVector();
2491          int * indices = regionSparse2->getIndices();
2492          int n = regionSparse2->getNumElements();
2493          memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double));
2494          double * array2 = save->denseVector();
2495          int * indices2 = save->getIndices();
2496          int n2 = save->getNumElements();
2497          assert (n == n2);
2498          if (save->packedMode()) {
2499               for (i = 0; i < n; i++) {
2500                    check[indices[i]] = array[i];
2501               }
2502               for (i = 0; i < n; i++) {
2503                    double value2 = array2[i];
2504                    assert (check[indices2[i]] == value2);
2505               }
2506          } else {
2507               int numberRows = coinFactorizationA_->numberRows();
2508               for (i = 0; i < numberRows; i++) {
2509                    double value1 = array[i];
2510                    double value2 = array2[i];
2511                    assert (value1 == value2);
2512               }
2513          }
2514          delete save;
2515          delete [] check;
2516          return returnCode;
2517#else
2518networkBasis_->updateColumn(regionSparse, regionSparse2, -1);
2519return 1;
2520#endif
2521     }
2522#endif
2523}
2524/* Updates one column (FTRAN) from region2
2525   Tries to do FT update
2526   number returned is negative if no room.
2527   Also updates region3
2528   region1 starts as zero and is zero at end */
2529int
2530ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
2531                                       CoinIndexedVector * regionSparse2,
2532                                       CoinIndexedVector * regionSparse3,
2533                                       bool noPermuteRegion3)
2534{
2535#ifdef CLP_DEBUG
2536     regionSparse1->checkClear();
2537#endif
2538     if (!numberRows())
2539          return 0;
2540     int returnCode = 0;
2541#ifndef SLIM_CLP
2542     if (!networkBasis_) {
2543#endif
2544#ifdef CLP_FACTORIZATION_INSTRUMENT
2545          factorization_instrument(-1);
2546#endif
2547          if (coinFactorizationA_) {
2548               coinFactorizationA_->setCollectStatistics(true);
2549               if (coinFactorizationA_->spaceForForrestTomlin()) {
2550                    assert (regionSparse2->packedMode());
2551                    assert (!regionSparse3->packedMode());
2552                    returnCode = coinFactorizationA_->updateTwoColumnsFT(regionSparse1,
2553                                 regionSparse2,
2554                                 regionSparse3,
2555                                 noPermuteRegion3);
2556               } else {
2557                    returnCode = coinFactorizationA_->updateColumnFT(regionSparse1,
2558                                 regionSparse2);
2559                    coinFactorizationA_->updateColumn(regionSparse1,
2560                                                      regionSparse3,
2561                                                      noPermuteRegion3);
2562               }
2563               coinFactorizationA_->setCollectStatistics(false);
2564          } else {
2565#if 0
2566               CoinSimpFactorization * fact =
2567                    dynamic_cast< CoinSimpFactorization*>(coinFactorizationB_);
2568               if (!fact) {
2569                    returnCode = coinFactorizationB_->updateColumnFT(regionSparse1,
2570                                 regionSparse2);
2571                    coinFactorizationB_->updateColumn(regionSparse1,
2572                                                      regionSparse3,
2573                                                      noPermuteRegion3);
2574               } else {
2575                    returnCode = fact->updateTwoColumnsFT(regionSparse1,
2576                                                          regionSparse2,
2577                                                          regionSparse3,
2578                                                          noPermuteRegion3);
2579               }
2580#else
2581#ifdef CLP_REUSE_ETAS
2582                    int tempInfo[2];
2583                    tempInfo[0] = model_->numberIterations();
2584                    tempInfo[1] = model_->sequenceIn();
2585                    coinFactorizationB_->setUsefulInformation(tempInfo, 3);
2586#endif
2587                    returnCode =
2588                      coinFactorizationB_->updateTwoColumnsFT(
2589                                                              regionSparse1,
2590                                                              regionSparse2,
2591                                                              regionSparse3,
2592                                                              noPermuteRegion3);
2593#endif
2594          }
2595#ifdef CLP_FACTORIZATION_INSTRUMENT
2596          factorization_instrument(9);
2597#endif
2598#ifdef PRINT_VECTOR
2599          printf("UpdateTwoFT\n");
2600          regionSparse2->print();
2601          regionSparse3->print();
2602#endif
2603          return returnCode;
2604#ifndef SLIM_CLP
2605     } else {
2606          returnCode = updateColumnFT(regionSparse1, regionSparse2);
2607          updateColumn(regionSparse1, regionSparse3, noPermuteRegion3);
2608     }
2609#endif
2610     return returnCode;
2611}
2612/* Updates one column (FTRAN) from region2
2613   number returned is negative if no room
2614   region1 starts as zero and is zero at end */
2615int
2616ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
2617          CoinIndexedVector * regionSparse2,
2618          bool noPermute) const
2619{
2620     if (!noPermute)
2621          regionSparse->checkClear();
2622     if (!coinFactorizationA_->numberRows())
2623          return 0;
2624     coinFactorizationA_->setCollectStatistics(false);
2625     int returnCode = coinFactorizationA_->updateColumn(regionSparse,
2626                      regionSparse2,
2627                      noPermute);
2628     return returnCode;
2629}
2630/* Updates one column (BTRAN) from region2
2631   region1 starts as zero and is zero at end */
2632int
2633ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
2634          CoinIndexedVector * regionSparse2) const
2635{
2636     if (!numberRows())
2637          return 0;
2638#ifndef SLIM_CLP
2639     if (!networkBasis_) {
2640#endif
2641#ifdef CLP_FACTORIZATION_INSTRUMENT
2642          factorization_instrument(-1);
2643#endif
2644          int returnCode;
2645
2646          if (coinFactorizationA_) {
2647               coinFactorizationA_->setCollectStatistics(true);
2648               returnCode =  coinFactorizationA_->updateColumnTranspose(regionSparse,
2649                             regionSparse2);
2650               coinFactorizationA_->setCollectStatistics(false);
2651          } else {
2652               returnCode = coinFactorizationB_->updateColumnTranspose(regionSparse,
2653                            regionSparse2);
2654          }
2655#ifdef CLP_FACTORIZATION_INSTRUMENT
2656          factorization_instrument(6);
2657#endif
2658#ifdef PRINT_VECTOR
2659          printf("UpdateTranspose\n");
2660          regionSparse2->print();
2661#endif
2662          return returnCode;
2663#ifndef SLIM_CLP
2664     } else {
2665#ifdef CHECK_NETWORK
2666          CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
2667          double * check = new double[coinFactorizationA_->numberRows()];
2668          int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
2669                           regionSparse2);
2670          networkBasis_->updateColumnTranspose(regionSparse, save);
2671          int i;
2672          double * array = regionSparse2->denseVector();
2673          int * indices = regionSparse2->getIndices();
2674          int n = regionSparse2->getNumElements();
2675          memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double));
2676          double * array2 = save->denseVector();
2677          int * indices2 = save->getIndices();
2678          int n2 = save->getNumElements();
2679          assert (n == n2);
2680          if (save->packedMode()) {
2681               for (i = 0; i < n; i++) {
2682                    check[indices[i]] = array[i];
2683               }
2684               for (i = 0; i < n; i++) {
2685                    double value2 = array2[i];
2686                    assert (check[indices2[i]] == value2);
2687               }
2688          } else {
2689               int numberRows = coinFactorizationA_->numberRows();
2690               for (i = 0; i < numberRows; i++) {
2691                    double value1 = array[i];
2692                    double value2 = array2[i];
2693                    assert (value1 == value2);
2694               }
2695          }
2696          delete save;
2697          delete [] check;
2698          return returnCode;
2699#else
2700return networkBasis_->updateColumnTranspose(regionSparse, regionSparse2);
2701#endif
2702     }
2703#endif
2704}
2705/* makes a row copy of L for speed and to allow very sparse problems */
2706void
2707ClpFactorization::goSparse()
2708{
2709#ifndef SLIM_CLP
2710     if (!networkBasis_) {
2711#endif
2712          if (coinFactorizationA_) {
2713#ifdef CLP_FACTORIZATION_INSTRUMENT
2714               factorization_instrument(-1);
2715#endif
2716               coinFactorizationA_->goSparse();
2717#ifdef CLP_FACTORIZATION_INSTRUMENT
2718               factorization_instrument(7);
2719#endif
2720          }
2721     }
2722}
2723// Cleans up i.e. gets rid of network basis
2724void
2725ClpFactorization::cleanUp()
2726{
2727#ifndef SLIM_CLP
2728     delete networkBasis_;
2729     networkBasis_ = NULL;
2730#endif
2731     if (coinFactorizationA_)
2732          coinFactorizationA_->resetStatistics();
2733}
2734/// Says whether to redo pivot order
2735bool
2736ClpFactorization::needToReorder() const
2737{
2738#ifdef CHECK_NETWORK
2739     return true;
2740#endif
2741#ifndef SLIM_CLP
2742     if (!networkBasis_)
2743#endif
2744          return true;
2745#ifndef SLIM_CLP
2746     else
2747          return false;
2748#endif
2749}
2750// Get weighted row list
2751void
2752ClpFactorization::getWeights(int * weights) const
2753{
2754#ifdef CLP_FACTORIZATION_INSTRUMENT
2755     factorization_instrument(-1);
2756#endif
2757#ifndef SLIM_CLP
2758     if (networkBasis_) {
2759          // Network - just unit
2760          int numberRows = coinFactorizationA_->numberRows();
2761          for (int i = 0; i < numberRows; i++)
2762               weights[i] = 1;
2763          return;
2764     }
2765#endif
2766     int * numberInRow = coinFactorizationA_->numberInRow();
2767     int * numberInColumn = coinFactorizationA_->numberInColumn();
2768     int * permuteBack = coinFactorizationA_->pivotColumnBack();
2769     int * indexRowU = coinFactorizationA_->indexRowU();
2770     const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
2771     const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
2772     int numberRows = coinFactorizationA_->numberRows();
2773     if (!startRowL || !coinFactorizationA_->numberInRow()) {
2774          int * temp = new int[numberRows];
2775          memset(temp, 0, numberRows * sizeof(int));
2776          int i;
2777          for (i = 0; i < numberRows; i++) {
2778               // one for pivot
2779               temp[i]++;
2780               CoinBigIndex j;
2781               for (j = startColumnU[i]; j < startColumnU[i] + numberInColumn[i]; j++) {
2782                    int iRow = indexRowU[j];
2783                    temp[iRow]++;
2784               }
2785          }
2786          CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
2787          int * indexRowL = coinFactorizationA_->indexRowL();
2788          int numberL = coinFactorizationA_->numberL();
2789          CoinBigIndex baseL = coinFactorizationA_->baseL();
2790          for (i = baseL; i < baseL + numberL; i++) {
2791               CoinBigIndex j;
2792               for (j = startColumnL[i]; j < startColumnL[i+1]; j++) {
2793                    int iRow = indexRowL[j];
2794                    temp[iRow]++;
2795               }
2796          }
2797          for (i = 0; i < numberRows; i++) {
2798               int number = temp[i];
2799               int iPermute = permuteBack[i];
2800               weights[iPermute] = number;
2801          }
2802          delete [] temp;
2803     } else {
2804          int i;
2805          for (i = 0; i < numberRows; i++) {
2806               int number = startRowL[i+1] - startRowL[i] + numberInRow[i] + 1;
2807               //number = startRowL[i+1]-startRowL[i]+1;
2808               //number = numberInRow[i]+1;
2809               int iPermute = permuteBack[i];
2810               weights[iPermute] = number;
2811          }
2812     }
2813#ifdef CLP_FACTORIZATION_INSTRUMENT
2814     factorization_instrument(8);
2815#endif
2816}
2817// Set tolerances to safer of existing and given
2818void
2819ClpFactorization::saferTolerances (  double zeroValue,
2820                                     double pivotValue)
2821{
2822     // better to have small tolerance even if slower
2823     zeroTolerance(CoinMin(zeroTolerance(), zeroValue));
2824     // better to have large tolerance even if slower
2825     double newValue;
2826     if (pivotValue > 0.0)
2827          newValue = pivotValue;
2828     else
2829          newValue = -pivotTolerance() * pivotValue;
2830     pivotTolerance(CoinMin(CoinMax(pivotTolerance(), newValue), 0.999));
2831}
2832// Sets factorization
2833void
2834ClpFactorization::setFactorization(ClpFactorization & rhs)
2835{
2836     ClpFactorization::operator=(rhs);
2837}
2838#endif
Note: See TracBrowser for help on using the repository browser.