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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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