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

Last change on this file since 2235 was 2235, checked in by forrest, 4 years ago

stuff for vector matrix

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