source: trunk/Clp/src/CoinAbcBaseFactorization5.cpp @ 2470

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

formatting

  • Property svn:keywords set to Id
File size: 27.7 KB
Line 
1/* $Id: CoinAbcBaseFactorization5.cpp 2385 2019-01-06 19:43:06Z stefan $ */
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#ifdef ABC_JUST_ONE_FACTORIZATION
6#include "CoinAbcCommonFactorization.hpp"
7#define CoinAbcTypeFactorization CoinAbcBaseFactorization
8#define ABC_SMALL -1
9#include "CoinAbcBaseFactorization.hpp"
10#endif
11#ifdef CoinAbcTypeFactorization
12
13#include "CoinIndexedVector.hpp"
14#include "CoinHelperFunctions.hpp"
15#include "CoinAbcHelperFunctions.hpp"
16#if ABC_NORMAL_DEBUG > 0
17// for conflicts
18extern int cilk_conflict;
19#endif
20#define UNROLL 1
21#define INLINE_IT
22//#define INLINE_IT2
23inline void scatterUpdateInline(CoinSimplexInt number,
24  CoinFactorizationDouble pivotValue,
25  const CoinFactorizationDouble *COIN_RESTRICT thisElement,
26  const CoinSimplexInt *COIN_RESTRICT thisIndex,
27  CoinFactorizationDouble *COIN_RESTRICT region)
28{
29#if UNROLL == 0
30  for (CoinBigIndex j = number - 1; j >= 0; j--) {
31    CoinSimplexInt iRow = thisIndex[j];
32    CoinFactorizationDouble regionValue = region[iRow];
33    CoinFactorizationDouble value = thisElement[j];
34    assert(value);
35    region[iRow] = regionValue - value * pivotValue;
36  }
37#elif UNROLL == 1
38  if ((number & 1) != 0) {
39    number--;
40    CoinSimplexInt iRow = thisIndex[number];
41    CoinFactorizationDouble regionValue = region[iRow];
42    CoinFactorizationDouble value = thisElement[number];
43    region[iRow] = regionValue - value * pivotValue;
44  }
45  for (CoinBigIndex j = number - 1; j >= 0; j -= 2) {
46    CoinSimplexInt iRow0 = thisIndex[j];
47    CoinSimplexInt iRow1 = thisIndex[j - 1];
48    CoinFactorizationDouble regionValue0 = region[iRow0];
49    CoinFactorizationDouble regionValue1 = region[iRow1];
50    region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
51    region[iRow1] = regionValue1 - thisElement[j - 1] * pivotValue;
52  }
53#elif UNROLL == 2
54  CoinSimplexInt iRow0;
55  CoinSimplexInt iRow1;
56  CoinFactorizationDouble regionValue0;
57  CoinFactorizationDouble regionValue1;
58  switch (static_cast< unsigned int >(number)) {
59  case 0:
60    break;
61  case 1:
62    iRow0 = thisIndex[0];
63    regionValue0 = region[iRow0];
64    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
65    break;
66  case 2:
67    iRow0 = thisIndex[0];
68    iRow1 = thisIndex[1];
69    regionValue0 = region[iRow0];
70    regionValue1 = region[iRow1];
71    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
72    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
73    break;
74  case 3:
75    iRow0 = thisIndex[0];
76    iRow1 = thisIndex[1];
77    regionValue0 = region[iRow0];
78    regionValue1 = region[iRow1];
79    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
80    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
81    iRow0 = thisIndex[2];
82    regionValue0 = region[iRow0];
83    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
84    break;
85  case 4:
86    iRow0 = thisIndex[0];
87    iRow1 = thisIndex[1];
88    regionValue0 = region[iRow0];
89    regionValue1 = region[iRow1];
90    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
91    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
92    iRow0 = thisIndex[2];
93    iRow1 = thisIndex[3];
94    regionValue0 = region[iRow0];
95    regionValue1 = region[iRow1];
96    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
97    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
98    break;
99  case 5:
100    iRow0 = thisIndex[0];
101    iRow1 = thisIndex[1];
102    regionValue0 = region[iRow0];
103    regionValue1 = region[iRow1];
104    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
105    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
106    iRow0 = thisIndex[2];
107    iRow1 = thisIndex[3];
108    regionValue0 = region[iRow0];
109    regionValue1 = region[iRow1];
110    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
111    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
112    iRow0 = thisIndex[4];
113    regionValue0 = region[iRow0];
114    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
115    break;
116  case 6:
117    iRow0 = thisIndex[0];
118    iRow1 = thisIndex[1];
119    regionValue0 = region[iRow0];
120    regionValue1 = region[iRow1];
121    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
122    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
123    iRow0 = thisIndex[2];
124    iRow1 = thisIndex[3];
125    regionValue0 = region[iRow0];
126    regionValue1 = region[iRow1];
127    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
128    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
129    iRow0 = thisIndex[4];
130    iRow1 = thisIndex[5];
131    regionValue0 = region[iRow0];
132    regionValue1 = region[iRow1];
133    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
134    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
135    break;
136  case 7:
137    iRow0 = thisIndex[0];
138    iRow1 = thisIndex[1];
139    regionValue0 = region[iRow0];
140    regionValue1 = region[iRow1];
141    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
142    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
143    iRow0 = thisIndex[2];
144    iRow1 = thisIndex[3];
145    regionValue0 = region[iRow0];
146    regionValue1 = region[iRow1];
147    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
148    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
149    iRow0 = thisIndex[4];
150    iRow1 = thisIndex[5];
151    regionValue0 = region[iRow0];
152    regionValue1 = region[iRow1];
153    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
154    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
155    iRow0 = thisIndex[6];
156    regionValue0 = region[iRow0];
157    region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
158    break;
159  case 8:
160    iRow0 = thisIndex[0];
161    iRow1 = thisIndex[1];
162    regionValue0 = region[iRow0];
163    regionValue1 = region[iRow1];
164    region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
165    region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
166    iRow0 = thisIndex[2];
167    iRow1 = thisIndex[3];
168    regionValue0 = region[iRow0];
169    regionValue1 = region[iRow1];
170    region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
171    region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
172    iRow0 = thisIndex[4];
173    iRow1 = thisIndex[5];
174    regionValue0 = region[iRow0];
175    regionValue1 = region[iRow1];
176    region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
177    region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
178    iRow0 = thisIndex[6];
179    iRow1 = thisIndex[7];
180    regionValue0 = region[iRow0];
181    regionValue1 = region[iRow1];
182    region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
183    region[iRow1] = regionValue1 - thisElement[7] * pivotValue;
184    break;
185  default:
186    if ((number & 1) != 0) {
187      number--;
188      CoinSimplexInt iRow = thisIndex[number];
189      CoinFactorizationDouble regionValue = region[iRow];
190      CoinFactorizationDouble value = thisElement[number];
191      region[iRow] = regionValue - value * pivotValue;
192    }
193    for (CoinBigIndex j = number - 1; j >= 0; j -= 2) {
194      CoinSimplexInt iRow0 = thisIndex[j];
195      CoinSimplexInt iRow1 = thisIndex[j - 1];
196      CoinFactorizationDouble regionValue0 = region[iRow0];
197      CoinFactorizationDouble regionValue1 = region[iRow1];
198      region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
199      region[iRow1] = regionValue1 - thisElement[j - 1] * pivotValue;
200    }
201    break;
202  }
203#endif
204}
205inline CoinFactorizationDouble gatherUpdate(CoinSimplexInt number,
206  const CoinFactorizationDouble *COIN_RESTRICT thisElement,
207  const CoinSimplexInt *COIN_RESTRICT thisIndex,
208  CoinFactorizationDouble *COIN_RESTRICT region)
209{
210  CoinFactorizationDouble pivotValue = 0.0;
211  for (CoinBigIndex j = 0; j < number; j++) {
212    CoinFactorizationDouble value = thisElement[j];
213    CoinSimplexInt jRow = thisIndex[j];
214    value *= region[jRow];
215    pivotValue -= value;
216  }
217  return pivotValue;
218}
219#undef INLINE_IT
220//  updateColumnR.  Updates part of column (FTRANR)
221void CoinAbcTypeFactorization::updateColumnR(CoinIndexedVector *regionSparse
222#if ABC_SMALL < 2
223  ,
224  CoinAbcStatistics &statistics
225#endif
226#if ABC_PARALLEL
227  ,
228  int whichSparse
229#endif
230  ) const
231{
232
233  if (numberR_) {
234    CoinSimplexDouble *COIN_RESTRICT region = regionSparse->denseVector();
235#if ABC_SMALL < 3
236    CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
237    CoinSimplexInt numberNonZero = regionSparse->getNumElements();
238#endif
239
240    const CoinBigIndex *COIN_RESTRICT startColumn = startColumnRAddress_ - numberRows_;
241    const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowRAddress_;
242    CoinSimplexInt *COIN_RESTRICT indexRowR2 = indexRowRAddress_ + lengthAreaR_;
243    if (gotRCopy())
244      indexRowR2 += lengthAreaR_;
245    const CoinFactorizationDouble *COIN_RESTRICT element = elementRAddress_;
246    const CoinSimplexInt *COIN_RESTRICT permute = permuteAddress_;
247    instrument_do("CoinAbcFactorizationUpdateR", sizeR);
248#if ABC_SMALL < 2
249    // Size of R
250    CoinSimplexDouble sizeR = startColumnRAddress_[numberR_];
251    // Work out very dubious idea of what would be fastest
252    // Average
253    CoinSimplexDouble averageR = sizeR / (static_cast< CoinSimplexDouble >(numberRowsExtra_));
254    // weights (relative to actual work)
255    CoinSimplexDouble setMark = 0.1; // setting mark
256    CoinSimplexDouble test1 = 1.0; // starting ftran (without testPivot)
257    CoinSimplexDouble testPivot = 2.0; // Seeing if zero etc
258    CoinSimplexDouble startDot = 2.0; // For starting dot product version
259    // For final scan
260    CoinSimplexDouble final = numberNonZero * 1.0;
261    CoinSimplexDouble methodTime0;
262    // For first type
263    methodTime0 = numberPivots_ * (testPivot + ((static_cast< double >(numberNonZero)) / (static_cast< double >(numberRows_)) * averageR));
264    methodTime0 += numberNonZero * (test1 + averageR);
265    methodTime0 += (numberNonZero + numberPivots_) * setMark;
266    // third
267    CoinSimplexDouble methodTime2 = sizeR + numberPivots_ * startDot + numberNonZero * final;
268    // switch off if necessary
269    CoinSimplexInt method = 0;
270    if (!gotRCopy() || !gotSparse() || methodTime2 < methodTime0) {
271      method = 2;
272    }
273    /*
274       IF we have row copy of R then always store in old way -
275       otherwise unpermuted way.
276    */
277#endif
278#if ABC_SMALL < 2
279    const CoinSimplexInt *numberInColumnPlus = numberInColumnPlusAddress_;
280#endif
281    const CoinSimplexInt *pivotRowBack = pivotColumnAddress_;
282#if ABC_SMALL < 2
283    switch (method) {
284    case 0: {
285      instrument_start("CoinAbcFactorizationUpdateRSparse1", numberRows_);
286      // mark known to be zero
287      CoinCheckZero *COIN_RESTRICT mark = markRowAddress_;
288#if ABC_PARALLEL
289      if (whichSparse) {
290        mark = reinterpret_cast< CoinCheckZero * >(reinterpret_cast< char * >(mark) + whichSparse * sizeSparseArray_);
291      }
292#endif
293      // mark all rows which will be permuted
294      for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++) {
295        CoinSimplexInt iRow = permute[i];
296        mark[iRow] = 1;
297      }
298      // we have another copy of R in R
299      const CoinFactorizationDouble *COIN_RESTRICT elementR = elementRAddress_ + lengthAreaR_;
300      const CoinSimplexInt *COIN_RESTRICT indexRowR = indexRowRAddress_ + lengthAreaR_;
301      const CoinBigIndex *COIN_RESTRICT startR = startColumnRAddress_ + maximumPivots_ + 1;
302      // For current list order does not matter as
303      // only affects end
304      CoinSimplexInt newNumber = 0;
305      for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
306        CoinSimplexInt iRow = regionIndex[i];
307        CoinFactorizationDouble pivotValue = region[iRow];
308        assert(region[iRow]);
309        if (!mark[iRow]) {
310          regionIndex[newNumber++] = iRow;
311        }
312        CoinSimplexInt kRow = permute[iRow];
313        CoinSimplexInt number = numberInColumnPlus[kRow];
314        instrument_add(number);
315        if (TEST_INT_NONZERO(number)) {
316          CoinBigIndex start = startR[kRow];
317#ifndef INLINE_IT
318          CoinBigIndex end = start + number;
319          for (CoinBigIndex j = start; j < end; j++) {
320            CoinFactorizationDouble value = elementR[j];
321            CoinSimplexInt jRow = indexRowR[j];
322            region[jRow] -= pivotValue * value;
323          }
324#else
325          CoinAbcScatterUpdate(number, pivotValue, elementR + start, indexRowR + start, region);
326#endif
327        }
328      }
329      instrument_start("CoinAbcFactorizationUpdateRSparse2", numberRows_);
330      numberNonZero = newNumber;
331      for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++) {
332        //move using permute_ (stored in inverse fashion)
333        CoinSimplexInt iRow = permute[i];
334        CoinFactorizationDouble pivotValue = region[iRow] + region[i];
335        //zero out pre-permuted
336        region[iRow] = 0.0;
337        region[i] = 0.0;
338        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
339          int jRow = pivotRowBack[i];
340          assert(!region[jRow]);
341          region[jRow] = pivotValue;
342          CoinSimplexInt number = numberInColumnPlus[i];
343          instrument_add(number);
344          if (TEST_INT_NONZERO(number)) {
345            CoinBigIndex start = startR[i];
346#ifndef INLINE_IT
347            CoinBigIndex end = start + number;
348            for (CoinBigIndex j = start; j < end; j++) {
349              CoinFactorizationDouble value = elementR[j];
350              CoinSimplexInt jRow = indexRowR[j];
351              region[jRow] -= pivotValue * value;
352            }
353#else
354            CoinAbcScatterUpdate(number, pivotValue, elementR + start, indexRowR + start, region);
355#endif
356          }
357        }
358      }
359      instrument_end();
360      for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++) {
361        CoinSimplexInt iRow = permute[i];
362        assert(iRow < numberRows_);
363        if (mark[iRow]) {
364          mark[iRow] = 0;
365          CoinExponent expValue = ABC_EXPONENT(region[iRow]);
366          if (expValue) {
367            if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
368              regionIndex[numberNonZero++] = iRow;
369            } else {
370              region[iRow] = 0.0;
371            }
372          }
373        }
374      }
375    } break;
376    case 2:
377#endif
378    {
379      instrument_start("CoinAbcFactorizationUpdateRDense", numberRows_);
380      CoinBigIndex start = startColumn[numberRows_];
381      for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++) {
382        //move using permute_ (stored in inverse fashion)
383        CoinBigIndex end = startColumn[i + 1];
384        CoinSimplexInt iRow = pivotRowBack[i];
385        assert(iRow < numberRows_);
386        CoinFactorizationDouble pivotValue = region[iRow];
387        instrument_add(end - start);
388        if (TEST_INT_NONZERO(end - start)) {
389#ifndef INLINE_IT2
390          for (CoinBigIndex j = start; j < end; j++) {
391            CoinFactorizationDouble value = element[j];
392            CoinSimplexInt jRow = indexRow[j];
393            value *= region[jRow];
394            pivotValue -= value;
395          }
396#else
397          pivotValue += CoinAbcGatherUpdate(end - start, element + start, indexRow + start, region);
398#endif
399          start = end;
400        }
401#if ABC_SMALL < 3
402        if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
403          if (!region[iRow])
404            regionIndex[numberNonZero++] = iRow;
405          region[iRow] = pivotValue;
406        } else {
407          if (region[iRow])
408            region[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
409        }
410#else
411        region[iRow] = pivotValue;
412#endif
413      }
414      instrument_end();
415#if ABC_SMALL < 3
416      // pack down
417      CoinSimplexInt n = numberNonZero;
418      numberNonZero = 0;
419      for (CoinSimplexInt i = 0; i < n; i++) {
420        CoinSimplexInt indexValue = regionIndex[i];
421        CoinExponent expValue = ABC_EXPONENT(region[indexValue]);
422        if (expValue) {
423          if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
424            regionIndex[numberNonZero++] = indexValue; // needed? is region used
425          } else {
426            region[indexValue] = 0.0;
427          }
428        }
429      }
430#endif
431    }
432#if ABC_SMALL < 2
433    break;
434    }
435#endif
436#if ABC_SMALL < 3
437    //set counts
438    regionSparse->setNumElements(numberNonZero);
439#endif
440  }
441#if ABC_SMALL < 2
442  if (factorizationStatistics())
443    statistics.countAfterR_ += regionSparse->getNumElements();
444#endif
445}
446#ifdef EARLY_FACTORIZE
447#include "AbcSimplex.hpp"
448#include "AbcMatrix.hpp"
449// Returns -2 if can't, -1 if singular, -99 memory, 0 OK
450int CoinAbcTypeFactorization::factorize(AbcSimplex *model, CoinIndexedVector &stuff)
451{
452  AbcMatrix *matrix = model->abcMatrix();
453  setStatus(-99);
454  int *COIN_RESTRICT pivotTemp = stuff.getIndices();
455  int numberColumnBasic = stuff.getNumElements();
456  // compute how much in basis
457  int numberSlacks = numberRows_ - numberColumnBasic;
458
459  CoinBigIndex numberElements = numberSlacks + matrix->countBasis(pivotTemp + numberSlacks, numberColumnBasic);
460  // Not needed for dense
461  numberElements = 3 * numberRows_ + 3 * numberElements + 20000;
462  getAreas(numberRows_, numberSlacks + numberColumnBasic, numberElements,
463    2 * numberElements);
464  numberSlacks_ = numberSlacks;
465  // Fill in counts so we can skip part of preProcess
466  // This is NOT needed for dense but would be needed for later versions
467  CoinFactorizationDouble *COIN_RESTRICT elementU;
468  int *COIN_RESTRICT indexRowU;
469  CoinBigIndex *COIN_RESTRICT startColumnU;
470  int *COIN_RESTRICT numberInRow;
471  int *COIN_RESTRICT numberInColumn;
472  elementU = this->elements();
473  indexRowU = this->indices();
474  startColumnU = this->starts();
475  numberInRow = this->numberInRow();
476  numberInColumn = this->numberInColumn();
477  CoinZeroN(numberInRow, numberRows_ + 1);
478  CoinZeroN(numberInColumn, numberRows_);
479  for (int i = 0; i < numberSlacks_; i++) {
480    int iRow = pivotTemp[i];
481    indexRowU[i] = iRow;
482    startColumnU[i] = i;
483    elementU[i] = 1.0;
484    numberInRow[iRow] = 1;
485    numberInColumn[i] = 1;
486  }
487  startColumnU[numberSlacks_] = numberSlacks_;
488  matrix->fillBasis(pivotTemp + numberSlacks_,
489    numberColumnBasic,
490    indexRowU,
491    startColumnU + numberSlacks_,
492    numberInRow,
493    numberInColumn + numberSlacks_,
494    elementU);
495  numberElements = startColumnU[numberRows_ - 1]
496    + numberInColumn[numberRows_ - 1];
497  preProcess();
498  factor(model);
499  if (status() == 0) {
500    // Put sequence numbers in workArea
501    int *savePivots = reinterpret_cast< int * >(workAreaAddress_);
502    CoinAbcMemcpy(savePivots, pivotTemp, numberRows_);
503    postProcess(pivotTemp, savePivots);
504    return 0;
505  } else {
506    return -2;
507  }
508}
509// 0 success, -1 can't +1 accuracy problems
510int CoinAbcTypeFactorization::replaceColumns(const AbcSimplex *model,
511  CoinIndexedVector &stuff,
512  int firstPivot, int lastPivot,
513  bool cleanUp)
514{
515  // Could skip some if goes in and out (only if easy as then no alpha check)
516  // Put sequence numbers in workArea
517  int *savePivots = reinterpret_cast< int * >(workAreaAddress_);
518  const int *pivotIndices = stuff.getIndices();
519  CoinSimplexDouble *pivotValues = stuff.denseVector();
520  int savePivot = stuff.capacity();
521  savePivot--;
522  savePivot -= 2 * firstPivot;
523  int numberDo = lastPivot - firstPivot;
524  // Say clear
525  stuff.setNumElements(0);
526  bool badUpdates = false;
527  for (int iPivot = 0; iPivot < numberDo; iPivot++) {
528    int sequenceOut = pivotIndices[--savePivot];
529    CoinSimplexDouble alpha = pivotValues[savePivot];
530    int sequenceIn = pivotIndices[--savePivot];
531    CoinSimplexDouble btranAlpha = pivotValues[savePivot];
532    int pivotRow;
533    for (pivotRow = 0; pivotRow < numberRows_; pivotRow++) {
534      if (sequenceOut == savePivots[pivotRow])
535        break;
536    }
537    assert(pivotRow < numberRows_);
538    savePivots[pivotRow] = sequenceIn;
539    model->unpack(stuff, sequenceIn);
540    updateColumnFTPart1(stuff);
541    stuff.clear();
542    //checkReplacePart1a(&stuff,pivotRow);
543    CoinSimplexDouble ftAlpha = checkReplacePart1(&stuff, pivotRow);
544    // may need better check
545    if (checkReplacePart2(pivotRow, btranAlpha, alpha, ftAlpha) > 1) {
546      badUpdates = true;
547      printf("Inaccuracy ? btranAlpha %g ftranAlpha %g ftAlpha %g\n",
548        btranAlpha, alpha, ftAlpha);
549      break;
550    }
551    replaceColumnPart3(model, &stuff, NULL, pivotRow, ftAlpha);
552  }
553  int flag;
554  if (!badUpdates) {
555    flag = 1;
556    if (cleanUp) {
557      CoinAbcMemcpy(model->pivotVariable(), savePivots, numberRows_);
558    }
559  } else {
560    flag = -1;
561  }
562  stuff.setNumElements(flag);
563  return flag > 0 ? 0 : flag;
564}
565#endif
566#ifdef ABC_ORDERED_FACTORIZATION
567// Permute in for Ftran
568void CoinAbcTypeFactorization::permuteInForFtran(CoinIndexedVector &regionSparse,
569  bool full) const
570{
571  int numberNonZero = regionSparse.getNumElements();
572  const int *COIN_RESTRICT permuteIn = permuteAddress_ + maximumRowsExtra_ + 1;
573  double *COIN_RESTRICT region = regionSparse.denseVector();
574  double *COIN_RESTRICT tempRegion = region + numberRows_;
575  int *COIN_RESTRICT index = regionSparse.getIndices();
576  if ((numberNonZero << 1) > numberRows_ || full) {
577    CoinAbcMemcpy(tempRegion, region, numberRows_);
578    CoinAbcMemset0(region, numberRows_);
579    numberNonZero = 0;
580    for (int i = 0; i < numberRows_; i++) {
581      double value = tempRegion[i];
582      if (value) {
583        tempRegion[i] = 0.0;
584        int iRow = permuteIn[i];
585        region[iRow] = value;
586        index[numberNonZero++] = iRow;
587      }
588    }
589    regionSparse.setNumElements(numberNonZero);
590  } else {
591    for (int i = 0; i < numberNonZero; i++) {
592      int iRow = index[i];
593      double value = region[iRow];
594      region[iRow] = 0.0;
595      iRow = permuteIn[iRow];
596      tempRegion[iRow] = value;
597      index[i] = iRow;
598    }
599    // and back
600    for (int i = 0; i < numberNonZero; i++) {
601      int iRow = index[i];
602      double value = tempRegion[iRow];
603      tempRegion[iRow] = 0.0;
604      region[iRow] = value;
605    }
606  }
607}
608// Permute in for Btran and multiply by pivot Region
609void CoinAbcTypeFactorization::permuteInForBtranAndMultiply(CoinIndexedVector &regionSparse,
610  bool full) const
611{
612  int numberNonZero = regionSparse.getNumElements();
613  double *COIN_RESTRICT region = regionSparse.denseVector();
614  double *COIN_RESTRICT tempRegion = region + numberRows_;
615  int *COIN_RESTRICT index = regionSparse.getIndices();
616  const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
617  if (full) {
618    numberNonZero = 0;
619    for (int iRow = 0; iRow < numberRows_; iRow++) {
620      double value = region[iRow];
621      if (value) {
622        region[iRow] = 0.0;
623        tempRegion[iRow] = value * pivotRegion[iRow];
624        index[numberNonZero++] = iRow;
625      }
626    }
627    regionSparse.setNumElements(numberNonZero);
628  } else {
629    for (int i = 0; i < numberNonZero; i++) {
630      int iRow = index[i];
631      double value = region[iRow];
632      region[iRow] = 0.0;
633      tempRegion[iRow] = value * pivotRegion[iRow];
634    }
635  }
636  regionSparse.setDenseVector(tempRegion);
637#if 0
638  if (numberNonZero)
639    printf("BtranIn vector %p region %p regionTemp %p\n",
640         &regionSparse,region,tempRegion);
641  else
642    printf("BtranIn no els! vector %p region %p regionTemp %p\n",
643         &regionSparse,region,tempRegion);
644#endif
645}
646// Permute out for Btran
647void CoinAbcTypeFactorization::permuteOutForBtran(CoinIndexedVector &regionSparse) const
648{
649  int numberNonZero = regionSparse.getNumElements();
650  const int *COIN_RESTRICT permuteIn = permuteAddress_ + maximumRowsExtra_ + 1;
651  const int *COIN_RESTRICT permuteOut = permuteIn + numberRows_;
652  CoinFactorizationDouble *COIN_RESTRICT tempRegion = denseVector(regionSparse);
653  CoinFactorizationDouble *COIN_RESTRICT region = tempRegion - numberRows_;
654  CoinSimplexInt *COIN_RESTRICT index = regionSparse.getIndices();
655  regionSparse.setDenseVector(region);
656  //printf("BtranOut vector %p region %p regionTemp %p\n",
657  //     &regionSparse,region,tempRegion);
658  for (int i = 0; i < numberNonZero; i++) {
659    int iRow = index[i];
660    double value = tempRegion[iRow];
661    tempRegion[iRow] = 0.0;
662    iRow = permuteOut[iRow];
663    region[iRow] = value;
664    index[i] = iRow;
665  }
666}
667#endif
668#if COIN_BIG_DOUBLE == 1
669// To a work array and associate vector
670void CoinAbcTypeFactorization::toLongArray(CoinIndexedVector *vector, int which) const
671{
672  assert(which >= 0 && which < FACTOR_CPU);
673  assert(!associatedVector_[which]);
674  associatedVector_[which] = vector;
675  CoinSimplexDouble *COIN_RESTRICT region = vector->denseVector();
676  CoinSimplexInt *COIN_RESTRICT regionIndex = vector->getIndices();
677  CoinSimplexInt numberNonZero = vector->getNumElements();
678  long double *COIN_RESTRICT longRegion = longArray_[which].array();
679  assert(!vector->packedMode());
680  // could check all of longRegion for zero but this should trap most errors
681  for (int i = 0; i < numberNonZero; i++) {
682    int iRow = regionIndex[i];
683    assert(!longRegion[iRow]);
684    longRegion[iRow] = region[iRow];
685    region[iRow] = 0.0;
686  }
687}
688// From a work array and dis-associate vector
689void CoinAbcTypeFactorization::fromLongArray(CoinIndexedVector *vector) const
690{
691  int which;
692  for (which = 0; which < FACTOR_CPU; which++)
693    if (vector == associatedVector_[which])
694      break;
695  assert(which < FACTOR_CPU);
696  associatedVector_[which] = NULL;
697  CoinSimplexDouble *COIN_RESTRICT region = vector->denseVector();
698  CoinSimplexInt *COIN_RESTRICT regionIndex = vector->getIndices();
699  CoinSimplexInt numberNonZero = vector->getNumElements();
700  long double *COIN_RESTRICT longRegion = longArray_[which].array();
701  assert(!vector->packedMode());
702  // could check all of region for zero but this should trap most errors
703  for (int i = 0; i < numberNonZero; i++) {
704    int iRow = regionIndex[i];
705    assert(!region[iRow]);
706    region[iRow] = longRegion[iRow];
707    longRegion[iRow] = 0.0;
708  }
709}
710// From a work array and dis-associate vector
711void CoinAbcTypeFactorization::fromLongArray(int which) const
712{
713  assert(which < FACTOR_CPU);
714  CoinIndexedVector *vector = associatedVector_[which];
715  associatedVector_[which] = NULL;
716  CoinSimplexDouble *COIN_RESTRICT region = vector->denseVector();
717  CoinSimplexInt *COIN_RESTRICT regionIndex = vector->getIndices();
718  CoinSimplexInt numberNonZero = vector->getNumElements();
719  long double *COIN_RESTRICT longRegion = longArray_[which].array();
720  assert(!vector->packedMode());
721  // could check all of region for zero but this should trap most errors
722  for (int i = 0; i < numberNonZero; i++) {
723    int iRow = regionIndex[i];
724    assert(!region[iRow]);
725    region[iRow] = longRegion[iRow];
726    longRegion[iRow] = 0.0;
727  }
728}
729// Returns long double * associated with vector
730long double *
731CoinAbcTypeFactorization::denseVector(CoinIndexedVector *vector) const
732{
733  int which;
734  for (which = 0; which < FACTOR_CPU; which++)
735    if (vector == associatedVector_[which])
736      break;
737  assert(which < FACTOR_CPU);
738  return longArray_[which].array();
739}
740// Returns long double * associated with vector
741long double *
742CoinAbcTypeFactorization::denseVector(CoinIndexedVector &vector) const
743{
744  int which;
745  for (which = 0; which < FACTOR_CPU; which++)
746    if (&vector == associatedVector_[which])
747      break;
748  assert(which < FACTOR_CPU);
749  return longArray_[which].array();
750}
751// Returns long double * associated with vector
752const long double *
753CoinAbcTypeFactorization::denseVector(const CoinIndexedVector *vector) const
754{
755  int which;
756  for (which = 0; which < FACTOR_CPU; which++)
757    if (vector == associatedVector_[which])
758      break;
759  assert(which < FACTOR_CPU);
760  return longArray_[which].array();
761}
762void CoinAbcTypeFactorization::scan(CoinIndexedVector *vector) const
763{
764  int numberNonZero = 0;
765  int *COIN_RESTRICT index = vector->getIndices();
766  CoinFactorizationDouble *COIN_RESTRICT region = denseVector(vector);
767  for (int i = 0; i < numberRows_; i++) {
768    if (fabs(region[i]) > zeroTolerance_)
769      index[numberNonZero++] = i;
770    else
771      region[i] = 0.0;
772  }
773  vector->setNumElements(numberNonZero);
774}
775// Clear all hidden arrays
776void CoinAbcTypeFactorization::clearHiddenArrays()
777{
778  for (int which = 0; which < FACTOR_CPU; which++) {
779    if (associatedVector_[which]) {
780      CoinIndexedVector *vector = associatedVector_[which];
781      associatedVector_[which] = NULL;
782      CoinFactorizationDouble *COIN_RESTRICT region = longArray_[which].array();
783      CoinSimplexInt *COIN_RESTRICT regionIndex = vector->getIndices();
784      CoinSimplexInt numberNonZero = vector->getNumElements();
785    }
786  }
787}
788#endif
789#endif
790
791/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
792*/
Note: See TracBrowser for help on using the repository browser.