source: stable/1.16/Clp/src/ClpFactorization.cpp @ 2116

Last change on this file since 2116 was 2116, checked in by forrest, 5 years ago

fix uninitialized

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