source: trunk/Clp/src/AbcSimplexFactorization.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 3 months ago

formatting

  • Property svn:keywords set to Id
File size: 32.8 KB
Line 
1/* $Id: AbcSimplexFactorization.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5#define USE_DENSE_FAC -1
6#define USE_SMALL_FAC 200
7#define USE_LONG_FAC 10000
8#include "CoinPragma.hpp"
9#include "AbcSimplexFactorization.hpp"
10#include "ClpFactorization.hpp"
11#include "ClpMessage.hpp"
12#include "CoinAbcCommon.hpp"
13#include "CoinHelperFunctions.hpp"
14#include "CoinIndexedVector.hpp"
15#include "AbcSimplex.hpp"
16#include "AbcSimplexDual.hpp"
17#include "AbcMatrix.hpp"
18#include "CoinAbcFactorization.hpp"
19#include "CoinFactorization.hpp"
20#ifdef ABC_JUST_ONE_FACTORIZATION
21#define CoinAbcFactorization CoinAbcBaseFactorization
22#define CoinAbcSmallFactorization CoinAbcBaseFactorization
23#define CoinAbcLongFactorization CoinAbcBaseFactorization
24#define CoinAbcOrderedFactorization CoinAbcBaseFactorization
25#endif
26#ifndef ABC_LONG_FACTORIZATION
27#undef CoinAbcLongFactorization
28#define CoinAbcLongFactorization CoinAbcOrderedFactorization
29#endif
30#ifdef ABC_TEMPORARY_FACTORIZATION
31#undef CoinAbcSmallFactorization
32#define CoinAbcSmallFactorization CoinAbcOrderedFactorization
33#endif
34#ifdef ABC_USE_COIN_FACTORIZATION
35#undef CoinAbcFactorization
36#undef CoinAbcSmallFactorization
37#undef CoinAbcLongFactorization
38#undef CoinAbcOrderedFactorization
39#define CoinAbcFactorization CoinFactorization
40#define CoinAbcSmallFactorization CoinFactorization
41#define CoinAbcLongFactorization CoinFactorization
42#define CoinAbcOrderedFactorization CoinFactorization
43#endif
44#ifndef ABC_USE_COIN_FACTORIZATION
45#undef CLP_FACTORIZATION_NEW_TIMING
46#else
47#ifndef CLP_FACTORIZATION_NEW_TIMING
48#define CLP_FACTORIZATION_NEW_TIMING 1
49#endif
50#endif
51//-------------------------------------------------------------------
52// Default Constructor
53//-------------------------------------------------------------------
54AbcSimplexFactorization::AbcSimplexFactorization(int /*numberRows*/)
55{
56  model_ = NULL;
57  coinAbcFactorization_ = new CoinAbcFactorization();
58  forceB_ = 0;
59  goDenseThreshold_ = USE_DENSE_FAC;
60  goSmallThreshold_ = USE_SMALL_FAC;
61  goLongThreshold_ = USE_LONG_FAC;
62  numberSlacks_ = 0;
63}
64
65//-------------------------------------------------------------------
66// Copy constructor
67//-------------------------------------------------------------------
68AbcSimplexFactorization::AbcSimplexFactorization(const AbcSimplexFactorization &rhs,
69  int denseIfSmaller)
70{
71  forceB_ = rhs.forceB_;
72  goDenseThreshold_ = rhs.goDenseThreshold_;
73  goSmallThreshold_ = rhs.goSmallThreshold_;
74  goLongThreshold_ = rhs.goLongThreshold_;
75  numberSlacks_ = rhs.numberSlacks_;
76  model_ = rhs.model_;
77#ifndef ABC_USE_COIN_FACTORIZATION
78  int goDense = 0;
79  if (denseIfSmaller > 0 && denseIfSmaller <= goDenseThreshold_) {
80    CoinAbcDenseFactorization *denseR = dynamic_cast< CoinAbcDenseFactorization * >(rhs.coinAbcFactorization_);
81    if (!denseR)
82      goDense = 1;
83  }
84  if (denseIfSmaller > 0 && !rhs.coinAbcFactorization_) {
85    if (denseIfSmaller <= goDenseThreshold_)
86      goDense = 1;
87    else if (denseIfSmaller <= goSmallThreshold_)
88      goDense = 2;
89    else if (denseIfSmaller >= goLongThreshold_)
90      goDense = 3;
91  } else if (denseIfSmaller < 0) {
92    if (-denseIfSmaller <= goDenseThreshold_)
93      goDense = 1;
94    else if (-denseIfSmaller <= goSmallThreshold_)
95      goDense = 2;
96    else if (-denseIfSmaller >= goLongThreshold_)
97      goDense = 3;
98  }
99  if (rhs.coinAbcFactorization_ && (denseIfSmaller >= 0 || !goDense))
100    coinAbcFactorization_ = rhs.coinAbcFactorization_->clone();
101  else
102    coinAbcFactorization_ = NULL;
103  if (goDense) {
104    delete coinAbcFactorization_;
105    if (goDense == 1)
106      coinAbcFactorization_ = new CoinAbcDenseFactorization();
107    else if (goDense == 2)
108      coinAbcFactorization_ = new CoinAbcSmallFactorization();
109    else if (goDense == 3)
110      coinAbcFactorization_ = new CoinAbcLongFactorization();
111    else
112      coinAbcFactorization_ = new CoinAbcFactorization();
113    assert(coinAbcFactorization_);
114    coinAbcFactorization_->maximumPivots(rhs.coinAbcFactorization_->maximumPivots());
115    coinAbcFactorization_->pivotTolerance(rhs.coinAbcFactorization_->pivotTolerance());
116    coinAbcFactorization_->zeroTolerance(rhs.coinAbcFactorization_->zeroTolerance());
117  }
118#else
119  coinAbcFactorization_ = new CoinFactorization(*rhs.coinAbcFactorization_);
120#endif
121}
122//-------------------------------------------------------------------
123// Destructor
124//-------------------------------------------------------------------
125AbcSimplexFactorization::~AbcSimplexFactorization()
126{
127  delete coinAbcFactorization_;
128}
129
130//----------------------------------------------------------------
131// Assignment operator
132//-------------------------------------------------------------------
133AbcSimplexFactorization &
134AbcSimplexFactorization::operator=(const AbcSimplexFactorization &rhs)
135{
136  if (this != &rhs) {
137    forceB_ = rhs.forceB_;
138    model_ = rhs.model_;
139    goDenseThreshold_ = rhs.goDenseThreshold_;
140    goSmallThreshold_ = rhs.goSmallThreshold_;
141    goLongThreshold_ = rhs.goLongThreshold_;
142    numberSlacks_ = rhs.numberSlacks_;
143
144    if (rhs.coinAbcFactorization_) {
145      delete coinAbcFactorization_;
146#ifndef ABC_USE_COIN_FACTORIZATION
147      coinAbcFactorization_ = rhs.coinAbcFactorization_->clone();
148#else
149      coinAbcFactorization_ = new CoinFactorization(*rhs.coinAbcFactorization_);
150#endif
151    }
152  }
153  return *this;
154}
155// Go over to dense code
156void AbcSimplexFactorization::goDenseOrSmall(int numberRows)
157{
158#ifndef ABC_USE_COIN_FACTORIZATION
159  if (!forceB_) {
160    delete coinAbcFactorization_;
161    if (numberRows <= goDenseThreshold_) {
162      coinAbcFactorization_ = new CoinAbcDenseFactorization();
163    } else if (numberRows <= goSmallThreshold_) {
164      coinAbcFactorization_ = new CoinAbcSmallFactorization();
165    } else if (numberRows >= goLongThreshold_) {
166      coinAbcFactorization_ = new CoinAbcLongFactorization();
167    } else {
168      coinAbcFactorization_ = new CoinAbcFactorization();
169    }
170  }
171#endif
172}
173// If nonzero force use of 1,dense 2,small 3,long
174void AbcSimplexFactorization::forceOtherFactorization(int which)
175{
176#ifndef ABC_USE_COIN_FACTORIZATION
177  delete coinAbcFactorization_;
178  forceB_ = 0;
179  coinAbcFactorization_ = NULL;
180  if (which > 0 && which < 6) {
181    forceB_ = which;
182    switch (which) {
183    case 1:
184      coinAbcFactorization_ = new CoinAbcDenseFactorization();
185      goDenseThreshold_ = COIN_INT_MAX;
186      break;
187    case 2:
188    case 4:
189      coinAbcFactorization_ = new CoinAbcSmallFactorization();
190      goSmallThreshold_ = COIN_INT_MAX;
191      break;
192    case 3:
193    case 5:
194      coinAbcFactorization_ = new CoinAbcLongFactorization();
195      goLongThreshold_ = 0;
196      break;
197    }
198  } else {
199    coinAbcFactorization_ = new CoinAbcFactorization();
200  }
201#endif
202}
203#ifdef CLP_FACTORIZATION_NEW_TIMING
204static bool readTwiddle = false;
205static double weightIncU = 1.0;
206static double weightR = 2.0;
207static double weightRest = 1.0;
208static double weightFactL = 30.0;
209static double weightFactDense = 0.1;
210static double weightNrows = 10.0;
211static double increaseNeeded = 1.1;
212static double constWeightIterate = 1.0;
213static double weightNrowsIterate = 3.0;
214bool AbcSimplexFactorization::timeToRefactorize() const
215{
216  bool reFactor = (coinAbcFactorization_->pivots() * 3 > coinAbcFactorization_->maximumPivots() * 2 && coinAbcFactorization_->numberElementsR() * 3 > (coinAbcFactorization_->numberElementsL() + coinAbcFactorization_->numberElementsU()) * 2 + 1000 && !coinAbcFactorization_->numberDense());
217  bool reFactor3 = false;
218  int numberPivots = coinAbcFactorization_->pivots();
219  //if (coinAbcFactorization_->pivots()<2)
220  if (numberPivots > lastNumberPivots_) {
221    if (!lastNumberPivots_) {
222      //lastR=0;
223      //lastU=endLengthU;
224      totalInR_ = 0.0;
225      totalInIncreasingU_ = 0.0;
226      shortestAverage_ = COIN_DBL_MAX;
227      if (!readTwiddle) {
228        readTwiddle = true;
229        char *environ = getenv("CLP_TWIDDLE");
230        if (environ) {
231          sscanf(environ, "%lg %lg %lg %lg %lg %lg %lg %lg %lg",
232            &weightIncU, &weightR, &weightRest, &weightFactL,
233            &weightFactDense, &weightNrows, &increaseNeeded,
234            &constWeightIterate, &weightNrowsIterate);
235        }
236        printf("weightIncU %g, weightR %g, weightRest %g, weightFactL %g, weightFactDense %g, weightNrows %g increaseNeeded %g constWeightIterate %g weightNrowsIterate %g\n",
237          weightIncU, weightR, weightRest, weightFactL,
238          weightFactDense, weightNrows, increaseNeeded,
239          constWeightIterate, weightNrowsIterate);
240      }
241    }
242    lastNumberPivots_ = numberPivots;
243    int numberDense = coinAbcFactorization_->numberDense();
244    double nnd = numberDense * numberDense;
245    int lengthL = coinAbcFactorization_->numberElementsL();
246    int lengthR = coinAbcFactorization_->numberElementsR();
247    int numberRows = coinAbcFactorization_->numberRows();
248    int lengthU = coinAbcFactorization_->numberElementsU() - (numberRows - numberDense);
249    totalInR_ += lengthR;
250    int effectiveU = lengthU - effectiveStartNumberU_;
251    totalInIncreasingU_ += effectiveU;
252    //lastR=lengthR;
253    //lastU=lengthU;
254    double rest = lengthL + 0.05 * nnd;
255    double constWeightFactor = weightFactL * lengthL + weightFactDense * nnd
256      + weightNrows * numberRows;
257    double constWeightIterateX = constWeightIterate * (lengthL + endLengthU_)
258      + weightNrowsIterate * numberRows;
259    double variableWeight = weightIncU * totalInIncreasingU_ + weightR * totalInR_ + weightRest * rest;
260    double average = constWeightIterateX + (constWeightFactor + variableWeight) / static_cast< double >(numberPivots);
261#if 0
262    if ((numberPivots%20)==0&&!ifPrint3)
263      printf("PIV %d nrow %d startU %d now %d L %d R %d dense %g average %g\n",
264             numberPivots,numberRows,effectiveStartNumberU_,
265             lengthU,lengthL,lengthR,nnd,average);
266#endif
267    shortestAverage_ = CoinMin(shortestAverage_, average);
268    if (average > increaseNeeded * shortestAverage_ && coinAbcFactorization_->pivots() > 30) {
269      //printf("PIVX %d nrow %d startU %d now %d L %d R %d dense %g average %g\n",
270      //     numberPivots,numberRows,effectiveStartNumberU_,
271      //     lengthU,lengthL,lengthR,nnd,average);
272      reFactor3 = true;
273    }
274  }
275  if (reFactor || reFactor3) {
276    reFactor = true;
277  }
278  return reFactor;
279}
280#if CLP_FACTORIZATION_NEW_TIMING > 1
281void AbcSimplexFactorization::statsRefactor(char when) const
282{
283  int numberPivots = coinAbcFactorization_->pivots();
284  int numberDense = coinAbcFactorization_->numberDense();
285  double nnd = numberDense * numberDense;
286  int lengthL = coinAbcFactorization_->numberElementsL();
287  int lengthR = coinAbcFactorization_->numberElementsR();
288  int numberRows = coinAbcFactorization_->numberRows();
289  int lengthU = coinAbcFactorization_->numberElementsU() - (numberRows - numberDense);
290  double rest = lengthL + 0.05 * nnd;
291  double constWeightFactor = weightFactL * lengthL + weightFactDense * nnd
292    + weightNrows * numberRows;
293  double constWeightIterateX = constWeightIterate * (lengthL + endLengthU_)
294    + weightNrowsIterate * numberRows;
295  double variableWeight = weightIncU * totalInIncreasingU_ + weightR * totalInR_ + weightRest * rest;
296  double average = constWeightIterateX + (constWeightFactor + variableWeight) / static_cast< double >(numberPivots);
297  printf("APIV%c %d nrow %d startU %d now %d L %d R %d dense %g average %g - shortest %g\n",
298    when, numberPivots, numberRows, effectiveStartNumberU_,
299    lengthU, lengthL, lengthR, nnd, average, shortestAverage_);
300}
301#endif
302#else
303bool AbcSimplexFactorization::timeToRefactorize() const
304{
305  return coinAbcFactorization_->pivots() > coinAbcFactorization_->numberRows() / 2.45 + 20;
306}
307#endif
308/* returns empty fake vector carved out of existing
309   later - maybe use associated arrays */
310static CoinIndexedVector *fakeVector(CoinIndexedVector *vector,
311  int fakeCapacity)
312{
313  int oldCapacity = vector->capacity();
314  CoinIndexedVector *newVector = new CoinIndexedVector();
315  newVector->setCapacity(fakeCapacity);
316  newVector->setDenseVector(vector->denseVector() + oldCapacity);
317  newVector->setIndexVector(vector->getIndices() + oldCapacity + ((oldCapacity + 3) >> 2));
318  vector->checkClean();
319  newVector->checkClear();
320  return newVector;
321}
322static void deleteFakeVector(CoinIndexedVector *vector,
323  CoinIndexedVector *fakeVector)
324{
325  int *index = vector->getIndices();
326  fakeVector->checkClear();
327  fakeVector->setCapacity(0);
328  fakeVector->setDenseVector(NULL);
329  fakeVector->setIndexVector(NULL);
330  delete fakeVector;
331  vector->checkClean();
332}
333// Synchronize stuff
334void AbcSimplexFactorization::synchronize(const ClpFactorization *otherFactorization, const AbcSimplex *model)
335{
336  goDenseThreshold_ = otherFactorization->goDenseThreshold();
337  goSmallThreshold_ = otherFactorization->goSmallThreshold();
338  goLongThreshold_ = otherFactorization->goOslThreshold();
339  //forceOtherFactorization(otherFactorization->typeOfFactorization());
340  goDenseOrSmall(model->numberRows());
341  maximumPivots(static_cast< int >(otherFactorization->maximumPivots() * 1.2));
342#ifdef ABC_USE_COIN_FACTORIZATION
343  // redo region sizes
344  int maximumRows = model->numberRows() + maximumPivots() + 1;
345  int currentCapacity = model->usefulArray(0)->capacity();
346  int newCapacity = currentCapacity + maximumRows + 3;
347  for (int i = 0; i < ABC_NUMBER_USEFUL; i++) {
348    CoinPartitionedVector *vector = model->usefulArray(i);
349    vector->reserve(newCapacity);
350    // zero
351    CoinZeroN(vector->getIndices(), newCapacity);
352    // but pretend old
353    //vector->setCapacity(currentCapacity);
354#if 0 //ndef NDEBUG
355    vector->checkClear();
356    CoinIndexedVector * newVector = fakeVector(vector,maximumRows);
357    deleteFakeVector(vector,newVector);
358#endif
359  }
360#endif
361}
362#ifdef CLP_FACTORIZATION_INSTRUMENT
363extern double externalTimeStart;
364extern double timeInFactorize;
365extern double timeInUpdate;
366extern double timeInUpdateTranspose;
367extern double timeInUpdateFT;
368extern double timeInUpdateTwoFT;
369extern double timeInReplace;
370extern int numberUpdate;
371extern int numberUpdateTranspose;
372extern int numberUpdateFT;
373extern int numberUpdateTwoFT;
374extern int numberReplace;
375extern int currentLengthR;
376extern int currentLengthU;
377extern int currentTakeoutU;
378#endif
379int AbcSimplexFactorization::factorize(AbcSimplex *model,
380  int solveType, bool valuesPass)
381{
382#ifdef CLP_FACTORIZATION_INSTRUMENT
383  externalTimeStart = CoinCpuTime();
384#endif
385  model_ = model;
386  AbcMatrix *matrix = model->abcMatrix();
387  int numberRows = model->numberRows();
388  if (!numberRows)
389    return 0;
390  bool anyChanged = false;
391  coinAbcFactorization_->setStatus(-99);
392#ifndef ABC_USE_COIN_FACTORIZATION
393  const
394#endif
395    int *COIN_RESTRICT pivotVariable
396    = model->pivotVariable();
397  //returns 0 -okay, -1 singular, -2 too many in basis */
398  // allow dense
399  int solveMode = coinAbcFactorization_->solveMode() & 1;
400  if (model->numberIterations() > model->baseIteration())
401    solveMode |= 9; // was +8 - this allows dense
402  else
403    solveMode = 1; // try dense
404  if (valuesPass)
405    solveMode += 4;
406  coinAbcFactorization_->setSolveMode(solveMode);
407  while (status() < -98) {
408
409    int i;
410    int numberBasic = 0;
411    // Move pivot variables across if they look good
412    int *COIN_RESTRICT pivotTemp = model->usefulArray(0)->getIndices();
413#ifndef NDEBUG
414    model_->checkArrays();
415#endif
416    assert(!model->usefulArray(0)->getNumElements());
417    // Seems to prefer things in order so quickest
418    // way is to go though like this
419    for (i = 0; i < numberRows; i++) {
420      if (pivotVariable[i] < numberRows)
421        pivotTemp[numberBasic++] = pivotVariable[i];
422    }
423    numberSlacks_ = numberBasic;
424    /* Put column basic variables into pivotVariable
425    */
426    for (i = 0; i < numberRows; i++) {
427      if (pivotVariable[i] >= numberRows)
428        pivotTemp[numberBasic++] = pivotVariable[i];
429    }
430    CoinBigIndex numberElements = numberSlacks_;
431
432    // compute how much in basis
433    int numberColumnBasic = numberBasic - numberSlacks_;
434
435    numberElements += matrix->countBasis(pivotTemp + numberSlacks_,
436      numberColumnBasic);
437    //printf("Basis has %d slacks - size %d\n",numberSlacks_,numberElements);
438    // Not needed for dense
439    numberElements = 3 * numberBasic + 3 * numberElements + 20000;
440#ifndef ABC_USE_COIN_FACTORIZATION
441    int numberIterations = model->numberIterations();
442    coinAbcFactorization_->setUsefulInformation(&numberIterations, 0);
443#else
444    coinAbcFactorization_->gutsOfDestructor();
445    coinAbcFactorization_->gutsOfInitialize(2);
446#endif
447    coinAbcFactorization_->getAreas(numberRows,
448      numberSlacks_ + numberColumnBasic, numberElements,
449      2 * numberElements);
450#if 0
451    if (!model->numberIterations())
452      printf("do I need destructor etc in getAreas?\n");
453#endif
454    // Fill in counts so we can skip part of preProcess
455    // This is NOT needed for dense but would be needed for later versions
456    CoinFactorizationDouble *COIN_RESTRICT elementU;
457    int *COIN_RESTRICT indexRowU;
458    CoinBigIndex *COIN_RESTRICT startColumnU;
459    int *COIN_RESTRICT numberInRow;
460    int *COIN_RESTRICT numberInColumn;
461#define slackValue 1.0
462#ifndef ABC_USE_COIN_FACTORIZATION
463    elementU = coinAbcFactorization_->elements();
464    indexRowU = coinAbcFactorization_->indices();
465    startColumnU = coinAbcFactorization_->starts();
466    numberInRow = coinAbcFactorization_->numberInRow();
467    numberInColumn = coinAbcFactorization_->numberInColumn();
468    coinAbcFactorization_->setNumberSlacks(numberSlacks_);
469    CoinZeroN(numberInRow, numberRows);
470    CoinZeroN(numberInColumn, numberRows);
471#else
472    elementU = coinAbcFactorization_->elementU();
473    indexRowU = coinAbcFactorization_->indexRowU();
474    startColumnU = coinAbcFactorization_->startColumnU();
475    numberInRow = coinAbcFactorization_->numberInRow();
476    numberInColumn = coinAbcFactorization_->numberInColumn();
477    CoinZeroN(numberInRow, coinAbcFactorization_->numberRows() + 1);
478    CoinZeroN(numberInColumn, coinAbcFactorization_->maximumColumnsExtra() + 1);
479#endif
480    for (i = 0; i < numberSlacks_; i++) {
481      int iRow = pivotTemp[i];
482      indexRowU[i] = iRow;
483      startColumnU[i] = i;
484      elementU[i] = slackValue;
485      numberInRow[iRow] = 1;
486      numberInColumn[i] = 1;
487    }
488    startColumnU[numberSlacks_] = numberSlacks_;
489    // can change for gub so redo
490    numberColumnBasic = numberRows - numberSlacks_;
491    matrix->fillBasis(pivotTemp + numberSlacks_,
492      numberColumnBasic,
493      indexRowU,
494      startColumnU + numberSlacks_,
495      numberInRow,
496      numberInColumn + numberSlacks_,
497      elementU);
498    numberElements = startColumnU[numberRows - 1]
499      + numberInColumn[numberRows - 1];
500#ifndef ABC_USE_COIN_FACTORIZATION
501    coinAbcFactorization_->preProcess();
502    coinAbcFactorization_->factor(model);
503#else
504    // recompute number basic
505    numberBasic = numberSlacks_ + numberColumnBasic;
506    if (numberBasic)
507      numberElements = startColumnU[numberBasic - 1]
508        + numberInColumn[numberBasic - 1];
509    else
510      numberElements = 0;
511    coinAbcFactorization_->setNumberElementsU(numberElements);
512#ifdef CLP_FACTORIZATION_NEW_TIMING
513    lastNumberPivots_ = 0;
514    effectiveStartNumberU_ = numberElements - numberRows;
515    //printf("%d slacks,%d in U at beginning\n",
516    //numberRowBasic,numberElements);
517#endif
518    //saveFactorization("dump.d");
519    if (coinAbcFactorization_->biasLU() >= 3 || coinAbcFactorization_->numberRows() != coinAbcFactorization_->numberColumns())
520      coinAbcFactorization_->preProcess(2);
521    else
522      coinAbcFactorization_->preProcess(3); // no row copy
523    coinAbcFactorization_->factor();
524#endif
525#if 0
526    if (model_->numberIterations()==23) {
527      CoinAbcFactorization * factor = dynamic_cast<CoinAbcFactorization *>(coinAbcFactorization_);
528      if (factor)
529        factor->show_self();
530    }
531#endif
532    if (coinAbcFactorization_->status() == -1 && (coinAbcFactorization_->solveMode() & 1) != 0) {
533      int solveMode = coinAbcFactorization_->solveMode();
534      solveMode--; // so bottom will be 0
535      coinAbcFactorization_->setSolveMode(solveMode);
536      coinAbcFactorization_->setStatus(-99);
537    } else if (coinAbcFactorization_->status() == -99) {
538      // get more memory
539      coinAbcFactorization_->areaFactor(coinAbcFactorization_->areaFactor() * 2.0);
540    }
541    if (coinAbcFactorization_->status() == -99)
542      continue;
543    // If we get here status is 0 or -1
544    if (coinAbcFactorization_->status() == 0 && numberBasic == numberRows) {
545#ifndef ABC_USE_COIN_FACTORIZATION
546      coinAbcFactorization_->postProcess(pivotTemp, model->pivotVariable());
547#else
548      const int *permuteBack = coinAbcFactorization_->permuteBack();
549      const int *back = coinAbcFactorization_->pivotColumnBack();
550      // Redo pivot order
551      for (i = 0; i < numberRows; i++) {
552        int k = pivotTemp[i];
553        // so rowIsBasic[k] would be permuteBack[back[i]]
554        int j = permuteBack[back[i]];
555        //assert (pivotVariable[j] == -1);
556        pivotVariable[j] = k;
557      }
558      // Set up permutation vector
559      // these arrays start off as copies of permute
560      // (and we could use permute_ instead of pivotColumn (not back though))
561      ClpDisjointCopyN(coinAbcFactorization_->permute(), numberRows, coinAbcFactorization_->pivotColumn());
562      ClpDisjointCopyN(coinAbcFactorization_->permuteBack(), numberRows, coinAbcFactorization_->pivotColumnBack());
563      // See if worth going sparse and when
564      coinAbcFactorization_->checkSparse();
565#ifdef CLP_FACTORIZATION_NEW_TIMING
566      endLengthU_ = coinAbcFactorization_->numberElements() - coinAbcFactorization_->numberDense() * coinAbcFactorization_->numberDense()
567        - coinAbcFactorization_->numberElementsL();
568#endif
569#endif
570      model_->moveToBasic();
571    } else if (solveType == 0 || solveType == 2 /*||solveType==1*/) {
572      // Change pivotTemp to be correct list
573      anyChanged = true;
574      coinAbcFactorization_->makeNonSingular(pivotTemp);
575      const double *COIN_RESTRICT lowerArray = model->lowerRegion();
576      const double *COIN_RESTRICT upperArray = model->upperRegion();
577      double *COIN_RESTRICT solution = model->solutionRegion();
578      //int * pivotVariable=model_->pivotVariable();
579      //int * fromExternal=model_->fromExternal();
580      int numberTotal = model_->numberTotal();
581      //can use external status_
582      unsigned char *COIN_RESTRICT statusArray = model_->statusArray();
583      CoinAbcMemset0(statusArray, numberTotal);
584      for (int iRow = 0; iRow < numberRows; iRow++) {
585        int iPivot = pivotVariable[iRow];
586        //if (iPivot!=pivotTemp[iRow])
587        //printf("row %d bad pivot %d good %d\n",iRow,iPivot,pivotTemp[iRow]);
588        statusArray[iPivot] = 1;
589      }
590      for (int iRow = 0; iRow < numberRows; iRow++) {
591        int iPivot = pivotTemp[iRow];
592        statusArray[iPivot] |= 2;
593      }
594      int jPivot = 0;
595      double largeValue = model->largeValue();
596      for (int iRow = 0; iRow < numberRows; iRow++) {
597        int iPivot = pivotVariable[iRow];
598        if (statusArray[iPivot] == 1) {
599          // clean
600          while (statusArray[jPivot] != 2) {
601            jPivot++;
602          }
603          statusArray[iPivot] = 0;
604          statusArray[jPivot] = 0;
605#ifndef NDEBUG
606          printf("On Row %d replacing %d by %d\n", iRow, iPivot, jPivot);
607#endif
608          AbcSimplex::Status thisStatus;
609          if (!valuesPass) {
610            double lower = lowerArray[iPivot];
611            double upper = upperArray[iPivot];
612            double value = solution[iPivot];
613            if (lower > -largeValue || upper < largeValue) {
614              if (lower != upper) {
615                if (fabs(value - lower) < fabs(value - upper)) {
616                  thisStatus = AbcSimplex::atLowerBound;
617                  solution[iPivot] = lower;
618                } else {
619                  thisStatus = AbcSimplex::atUpperBound;
620                  solution[iPivot] = upper;
621                }
622              } else {
623                thisStatus = AbcSimplex::isFixed;
624              }
625            } else {
626              thisStatus = AbcSimplex::isFree;
627            }
628          } else {
629            thisStatus = AbcSimplex::superBasic;
630          }
631          model_->setInternalStatus(iPivot, thisStatus);
632          model_->setInternalStatus(jPivot, AbcSimplex::basic);
633          // swap (solution will be wrong - but that doesn't matter as basic)
634          model_->swap(iRow, jPivot);
635        }
636      }
637      CoinAbcMemcpy(model_->pivotVariable(), pivotTemp, numberRows);
638#ifndef NDEBUG
639      model_->checkConsistentPivots();
640#endif
641      // signal repeat
642      coinAbcFactorization_->pivotTolerance(0.999);
643      coinAbcFactorization_->setStatus(-99);
644    }
645  }
646  //coinAbcFactorization_->setSolveMode(solveMode|1);
647  if (anyChanged && model->algorithm() < 0 && solveType > 0) {
648    double dummyCost;
649    static_cast< AbcSimplexDual * >(model)->changeBounds(3, dummyCost);
650  }
651  return coinAbcFactorization_->status();
652}
653#if 0
654/* Checks if can replace one Column to basis,
655   returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots
656   Fills in region for use later
657   partial update already in U */
658int
659AbcSimplexFactorization::checkReplace ( const AbcSimplex * model,
660                                        CoinIndexedVector * regionSparse,
661                                        int pivotRow,
662                                        double &pivotCheck,
663                                        double acceptablePivot)
664{
665  if (pivots()==maximumPivots())
666    return 5;
667  else
668    return coinAbcFactorization_->checkReplace(model,regionSparse,pivotRow,pivotCheck,acceptablePivot);
669}
670/* Replaces one Column in basis,
671   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
672   If skipBtranU is false will do btran part
673   partial update already in U */
674int
675AbcSimplexFactorization::replaceColumn ( const AbcSimplex * model,
676                                         CoinIndexedVector * regionSparse,
677                                         CoinIndexedVector * tableauColumn,
678                                         int pivotRow,
679                                         double pivotCheck ,
680                                         bool skipBtranU,
681                                         double acceptablePivot)
682{
683  bool tab = coinAbcFactorization_->wantsTableauColumn();
684  int tempInfo[1];
685  tempInfo[0] = model->numberIterations();
686  coinAbcFactorization_->setUsefulInformation(tempInfo, 1);
687#ifdef CLP_FACTORIZATION_NEW_TIMING
688  int nOld=0;
689  int nNew=0;
690  int seq;
691  const CoinPackedMatrix * matrix=model->matrix();
692  const int * columnLength = matrix->getVectorLengths();
693  seq=model->sequenceIn();
694  if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
695    if (seq<model->numberRows())
696      nNew=1;
697    else
698      nNew=columnLength[seq-model->numberRows()];
699  }
700  seq=model->sequenceOut();
701  if (seq>=0&&seq<model->numberColumns()+model->numberRows()) {
702    if (seq<model->numberRows())
703      nOld=1;
704    else
705      nOld=columnLength[seq-model->numberRows()];
706  }
707  effectiveStartNumberU_ += nNew-nOld;
708#endif
709  int returnCode =
710    coinAbcFactorization_->replaceColumn(tab ? tableauColumn : regionSparse,
711                                       pivotRow,
712                                       pivotCheck,
713                                       skipBtranU,
714                                       acceptablePivot);
715  return returnCode;
716}
717#endif
718/* Replaces one Column to basis,
719   partial update already in U */
720void AbcSimplexFactorization::replaceColumnPart3(const AbcSimplex *model,
721  CoinIndexedVector *regionSparse,
722  CoinIndexedVector *tableauColumn,
723  int pivotRow,
724#ifdef ABC_LONG_FACTORIZATION
725  long
726#endif
727  double alpha)
728{
729#ifndef ABC_USE_COIN_FACTORIZATION
730  bool tab = coinAbcFactorization_->wantsTableauColumn();
731  int tempInfo[1];
732  tempInfo[0] = model->numberIterations();
733  coinAbcFactorization_->setUsefulInformation(tempInfo, 1);
734  if (tab)
735    coinAbcFactorization_->replaceColumnPart3(model, NULL, tableauColumn,
736      pivotRow,
737      tableauColumn->denseVector()[pivotRow]);
738  else
739    coinAbcFactorization_->replaceColumnPart3(model, regionSparse, NULL,
740      pivotRow,
741      alpha);
742#else
743#ifdef CLP_FACTORIZATION_NEW_TIMING
744  int nOld = 0;
745  int nNew = 0;
746  int seq;
747  const CoinPackedMatrix *matrix = model->matrix();
748  const int *columnLength = matrix->getVectorLengths();
749  seq = model->sequenceIn();
750  if (seq >= 0 && seq < model->numberColumns() + model->numberRows()) {
751    if (seq < model->numberRows())
752      nNew = 1;
753    else
754      nNew = columnLength[seq - model->numberRows()];
755  }
756  seq = model->sequenceOut();
757  if (seq >= 0 && seq < model->numberColumns() + model->numberRows()) {
758    if (seq < model->numberRows())
759      nOld = 1;
760    else
761      nOld = columnLength[seq - model->numberRows()];
762  }
763  effectiveStartNumberU_ += nNew - nOld;
764#endif
765  coinAbcFactorization_->replaceColumnPart3(regionSparse,
766    pivotRow,
767    alpha);
768#endif
769}
770/* Replaces one Column to basis,
771   partial update in vector */
772void AbcSimplexFactorization::replaceColumnPart3(const AbcSimplex *model,
773  CoinIndexedVector *regionSparse,
774  CoinIndexedVector *tableauColumn,
775  CoinIndexedVector *partialUpdate,
776  int pivotRow,
777#ifdef ABC_LONG_FACTORIZATION
778  long
779#endif
780  double alpha)
781{
782#ifndef ABC_USE_COIN_FACTORIZATION
783  bool tab = coinAbcFactorization_->wantsTableauColumn();
784  int tempInfo[1];
785  tempInfo[0] = model->numberIterations();
786  coinAbcFactorization_->setUsefulInformation(tempInfo, 1);
787  if (tab)
788    coinAbcFactorization_->replaceColumnPart3(model, NULL, tableauColumn,
789      pivotRow,
790      tableauColumn->denseVector()[pivotRow]);
791  else
792    coinAbcFactorization_->replaceColumnPart3(model, regionSparse, NULL, partialUpdate,
793      pivotRow,
794      alpha);
795#else
796  coinAbcFactorization_->replaceColumnPart3(regionSparse, partialUpdate,
797    pivotRow,
798    alpha);
799#endif
800}
801#if 0
802/* Updates one column (FTRAN) from region2
803   number returned is negative if no room
804   region1 starts as zero and is zero at end */
805int
806AbcSimplexFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
807                                          CoinIndexedVector * regionSparse2)
808{
809  if (!numberRows())
810    return 0;
811  int returnCode;
812  returnCode = coinAbcFactorization_->updateColumnFT(regionSparse,
813                                                   regionSparse2);
814  return returnCode;
815}
816/* Updates one column (FTRAN) from region2
817   number returned is negative if no room
818   region1 starts as zero and is zero at end */
819int
820AbcSimplexFactorization::updateColumn ( CoinIndexedVector * regionSparse,
821                                        CoinIndexedVector * regionSparse2) const
822{
823  if (!numberRows())
824    return 0;
825  int returnCode;
826  returnCode = coinAbcFactorization_->updateColumn(regionSparse,
827                                                 regionSparse2,true);
828  return returnCode;
829}
830/* Updates one column (FTRAN) from region2
831   Tries to do FT update
832   number returned is negative if no room.
833   Also updates region3
834   region1 starts as zero and is zero at end */
835int
836AbcSimplexFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
837                                              CoinIndexedVector * regionSparse2,
838                                              CoinIndexedVector * regionSparse3)
839{
840  if (!numberRows())
841    return 0;
842  int returnCode = 0;
843  returnCode =
844    coinAbcFactorization_->updateTwoColumnsFT(
845                                            regionSparse1,
846                                            regionSparse2,
847                                            regionSparse3,
848                                            true);
849  return returnCode;
850}
851/* Updates one column (BTRAN) from region2
852   region1 starts as zero and is zero at end */
853int
854AbcSimplexFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
855                                                 CoinIndexedVector * regionSparse2) const
856{
857  if (!numberRows())
858    return 0;
859  int returnCode;
860  returnCode = coinAbcFactorization_->updateColumnTranspose(regionSparse,
861                                                          regionSparse2);
862  return returnCode;
863}
864/* Updates one column for dual steepest edge weights (FTRAN) */
865void
866AbcSimplexFactorization::updateWeights ( CoinIndexedVector & regionSparse) const
867{
868  // NOTE either switch off sparse or pass in a sparseArray_ so can go parallel
869  // may be best to use inner product approach
870  static double fraction[2]={0.0,0.0};
871  static int times=0;
872  times++;
873  fraction[0] += static_cast<double>(regionSparse.getNumElements())/
874    (static_cast<double>(model_->numberRows())+1.0);
875  updateColumn(regionSparse);
876  fraction[1] += static_cast<double>(regionSparse.getNumElements())/
877    (static_cast<double>(model_->numberRows())+1.0);
878  if ((times%1000)==0)
879    printf("Average density %g before then %g\n",
880           (100.0*fraction[0])/static_cast<double>(times),
881           (100.0*fraction[1])/static_cast<double>(times));
882}
883#endif
884/* makes a row copy of L for speed and to allow very sparse problems */
885void AbcSimplexFactorization::goSparse()
886{
887#ifdef ABC_USE_COIN_FACTORIZATION
888  // sparse methods
889  coinAbcFactorization_->sparseThreshold(0);
890  coinAbcFactorization_->goSparse();
891#else
892  abort();
893  coinAbcFactorization_->goSparse();
894#endif
895}
896// Set tolerances to safer of existing and given
897void AbcSimplexFactorization::saferTolerances(double zeroValue,
898  double pivotValue)
899{
900  double newValue1;
901  // better to have small tolerance even if slower
902  if (zeroValue > 0.0)
903    newValue1 = zeroValue;
904  else
905    newValue1 = -zeroTolerance() * zeroValue;
906  newValue1 = CoinMin(zeroTolerance(), newValue1);
907  if (newValue1 > 1.0e-15)
908    zeroTolerance(newValue1);
909  double newValue2;
910  // better to have large tolerance even if slower
911  if (pivotValue > 0.0)
912    newValue2 = pivotValue;
913  else
914    newValue2 = -pivotTolerance() * pivotValue;
915  newValue2 = CoinMin(CoinMax(pivotTolerance(), newValue2), 0.999);
916  if (newValue2 > pivotTolerance()) {
917    pivotTolerance(newValue2);
918    char line[100];
919    sprintf(line, "new zero tolerance %g new pivot tolerance %g",
920      zeroTolerance(), pivotTolerance());
921    model_->messageHandler()->message(CLP_GENERAL2, *model_->messagesPointer())
922      << line << CoinMessageEol;
923  }
924}
925// Sets factorization
926void AbcSimplexFactorization::setFactorization(AbcSimplexFactorization &rhs)
927{
928  AbcSimplexFactorization::operator=(rhs);
929}
930
931/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
932*/
Note: See TracBrowser for help on using the repository browser.