source: trunk/Clp/src/CoinAbcDenseFactorization.cpp @ 2449

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

formatting

  • Property svn:keywords set to Id
File size: 24.0 KB
Line 
1/* $Id: CoinAbcDenseFactorization.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2008, 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
6#include "CoinUtilsConfig.h"
7#include "CoinPragma.hpp"
8
9#include <cassert>
10#include <cstdio>
11
12#include "CoinAbcDenseFactorization.hpp"
13#include "CoinAbcCommonFactorization.hpp"
14#include "CoinIndexedVector.hpp"
15#include "CoinHelperFunctions.hpp"
16#include "CoinPackedMatrix.hpp"
17#include "CoinFinite.hpp"
18//:class CoinAbcDenseFactorization.  Deals with Factorization and Updates
19//  CoinAbcDenseFactorization.  Constructor
20CoinAbcDenseFactorization::CoinAbcDenseFactorization()
21  : CoinAbcAnyFactorization()
22{
23  gutsOfInitialize();
24}
25
26/// Copy constructor
27CoinAbcDenseFactorization::CoinAbcDenseFactorization(const CoinAbcDenseFactorization &other)
28  : CoinAbcAnyFactorization(other)
29{
30  gutsOfInitialize();
31  gutsOfCopy(other);
32}
33// Clone
34CoinAbcAnyFactorization *
35CoinAbcDenseFactorization::clone() const
36{
37  return new CoinAbcDenseFactorization(*this);
38}
39/// The real work of constructors etc
40void CoinAbcDenseFactorization::gutsOfDestructor()
41{
42  delete[] elements_;
43  delete[] pivotRow_;
44  delete[] workArea_;
45  elements_ = NULL;
46  pivotRow_ = NULL;
47  workArea_ = NULL;
48  numberRows_ = 0;
49  numberGoodU_ = 0;
50  status_ = -1;
51  maximumRows_ = 0;
52#if ABC_PARALLEL == 2
53  parallelMode_ = 0;
54#endif
55  maximumSpace_ = 0;
56  maximumRowsAdjusted_ = 0;
57  solveMode_ = 0;
58}
59void CoinAbcDenseFactorization::gutsOfInitialize()
60{
61  pivotTolerance_ = 1.0e-1;
62  minimumPivotTolerance_ = 1.0e-1;
63  zeroTolerance_ = 1.0e-13;
64  areaFactor_ = 1.0;
65  numberDense_ = 0;
66  maximumPivots_ = 200;
67  relaxCheck_ = 1.0;
68  numberRows_ = 0;
69  numberGoodU_ = 0;
70  status_ = -1;
71  numberSlacks_ = 0;
72  numberPivots_ = 0;
73  maximumRows_ = 0;
74#if ABC_PARALLEL == 2
75  parallelMode_ = 0;
76#endif
77  maximumSpace_ = 0;
78  maximumRowsAdjusted_ = 0;
79  elements_ = NULL;
80  pivotRow_ = NULL;
81  workArea_ = NULL;
82  solveMode_ = 0;
83}
84//  ~CoinAbcDenseFactorization.  Destructor
85CoinAbcDenseFactorization::~CoinAbcDenseFactorization()
86{
87  gutsOfDestructor();
88}
89//  =
90CoinAbcDenseFactorization &CoinAbcDenseFactorization::operator=(const CoinAbcDenseFactorization &other)
91{
92  if (this != &other) {
93    gutsOfDestructor();
94    gutsOfInitialize();
95    gutsOfCopy(other);
96  }
97  return *this;
98}
99#define WORK_MULT 2
100void CoinAbcDenseFactorization::gutsOfCopy(const CoinAbcDenseFactorization &other)
101{
102  pivotTolerance_ = other.pivotTolerance_;
103  minimumPivotTolerance_ = other.minimumPivotTolerance_;
104  zeroTolerance_ = other.zeroTolerance_;
105  areaFactor_ = other.areaFactor_;
106  numberDense_ = other.numberDense_;
107  relaxCheck_ = other.relaxCheck_;
108  numberRows_ = other.numberRows_;
109  maximumRows_ = other.maximumRows_;
110#if ABC_PARALLEL == 2
111  parallelMode_ = other.parallelMode_;
112#endif
113  maximumSpace_ = other.maximumSpace_;
114  maximumRowsAdjusted_ = other.maximumRowsAdjusted_;
115  solveMode_ = other.solveMode_;
116  numberGoodU_ = other.numberGoodU_;
117  maximumPivots_ = other.maximumPivots_;
118  numberPivots_ = other.numberPivots_;
119  factorElements_ = other.factorElements_;
120  status_ = other.status_;
121  numberSlacks_ = other.numberSlacks_;
122  if (other.pivotRow_) {
123    pivotRow_ = new int[2 * maximumRowsAdjusted_ + maximumPivots_];
124    CoinMemcpyN(other.pivotRow_, (2 * maximumRowsAdjusted_ + numberPivots_), pivotRow_);
125    elements_ = new CoinFactorizationDouble[maximumSpace_];
126    CoinMemcpyN(other.elements_, (maximumRowsAdjusted_ + numberPivots_) * maximumRowsAdjusted_, elements_);
127    workArea_ = new CoinFactorizationDouble[maximumRowsAdjusted_ * WORK_MULT];
128    CoinZeroN(workArea_, maximumRowsAdjusted_ * WORK_MULT);
129  } else {
130    elements_ = NULL;
131    pivotRow_ = NULL;
132    workArea_ = NULL;
133  }
134}
135
136//  getAreas.  Gets space for a factorization
137//called by constructors
138void CoinAbcDenseFactorization::getAreas(int numberOfRows,
139  int /*numberOfColumns*/,
140  CoinBigIndex,
141  CoinBigIndex)
142{
143
144  numberRows_ = numberOfRows;
145  numberDense_ = numberRows_;
146  if ((numberRows_ & (BLOCKING8 - 1)) != 0)
147    numberDense_ += (BLOCKING8 - (numberRows_ & (BLOCKING8 - 1)));
148  CoinBigIndex size = numberDense_ * (2 * numberDense_ + CoinMax(maximumPivots_, (numberDense_ + 1) >> 1));
149  if (size > maximumSpace_) {
150    delete[] elements_;
151    elements_ = new CoinFactorizationDouble[size];
152    maximumSpace_ = size;
153  }
154  if (numberRows_ > maximumRows_) {
155    maximumRows_ = numberRows_;
156    maximumRowsAdjusted_ = maximumRows_;
157    if ((maximumRows_ & (BLOCKING8 - 1)) != 0)
158      maximumRowsAdjusted_ += (BLOCKING8 - (maximumRows_ & (BLOCKING8 - 1)));
159    delete[] pivotRow_;
160    delete[] workArea_;
161    pivotRow_ = new int[2 * maximumRowsAdjusted_ + maximumPivots_];
162    workArea_ = new CoinFactorizationDouble[maximumRowsAdjusted_ * WORK_MULT];
163  }
164}
165
166//  preProcess.
167void CoinAbcDenseFactorization::preProcess()
168{
169  CoinBigIndex put = numberDense_ * numberDense_;
170  CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
171  CoinAbcMemset0(area, put);
172  int *indexRow = reinterpret_cast< int * >(elements_ + numberRows_ * numberRows_);
173  CoinBigIndex *starts = reinterpret_cast< CoinBigIndex * >(pivotRow_);
174  CoinFactorizationDouble *COIN_RESTRICT column = area;
175  //if (solveMode_==10)
176  //solveMode_=1;
177  if ((solveMode_ % 10) != 0) {
178    for (int i = 0; i < numberRows_; i++) {
179      for (CoinBigIndex j = starts[i]; j < starts[i + 1]; j++) {
180        int iRow = indexRow[j];
181        int iBlock = iRow >> 3;
182        iRow = iRow & (BLOCKING8 - 1);
183        column[iRow + BLOCKING8X8 * iBlock] = elements_[j];
184      }
185      if ((i & (BLOCKING8 - 1)) != (BLOCKING8 - 1))
186        column += BLOCKING8;
187      else
188        column += numberDense_ * BLOCKING8 - (BLOCKING8 - 1) * BLOCKING8;
189    }
190    for (int i = numberRows_; i < numberDense_; i++) {
191      int iRow = i;
192      int iBlock = iRow >> 3;
193      iRow = iRow & (BLOCKING8 - 1);
194      column[iRow + BLOCKING8X8 * iBlock] = 1.0;
195      if ((i & (BLOCKING8 - 1)) != (BLOCKING8 - 1))
196        column += BLOCKING8;
197      //else
198      //column += numberDense_*BLOCKING8-(BLOCKING8-1)*BLOCKING8;
199    }
200  } else {
201    // first or bad
202    for (int i = 0; i < numberRows_; i++) {
203      for (CoinBigIndex j = starts[i]; j < starts[i + 1]; j++) {
204        int iRow = indexRow[j];
205        column[iRow] = elements_[j];
206      }
207      column += numberDense_;
208    }
209    for (int i = numberRows_; i < numberDense_; i++) {
210      column[i] = 1.0;
211      column += numberDense_;
212    }
213  }
214}
215
216//Does factorization
217int CoinAbcDenseFactorization::factor(AbcSimplex * /*model*/)
218{
219  numberPivots_ = 0;
220  // ? take out
221  //printf("Debug non lapack dense factor\n");
222  //solveMode_&=~1;
223  CoinBigIndex put = numberDense_ * numberDense_;
224  CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
225  if ((solveMode_ % 10) != 0) {
226    // save last start
227    CoinBigIndex lastStart = pivotRow_[numberRows_];
228    status_ = CoinAbcDgetrf(numberDense_, numberDense_, area, numberDense_, pivotRow_ + numberDense_
229#if ABC_PARALLEL == 2
230      ,
231      parallelMode_
232#endif
233    );
234    // need to check size of pivots
235    if (!status_) {
236      // OK
237      solveMode_ = 1 + 10 * (solveMode_ / 10);
238      numberGoodU_ = numberRows_;
239      CoinZeroN(workArea_, 2 * numberRows_);
240#if 0 //ndef NDEBUG
241      const CoinFactorizationDouble * column = area;
242      double smallest=COIN_DBL_MAX;
243      for (int i=0;i<numberRows_;i++) {
244        if (fabs(column[i])<smallest)
245          smallest = fabs(column[i]);
246        column += numberRows_;
247      }
248      if (smallest<1.0e-8)
249        printf("small el %g\n",smallest);
250#endif
251      return 0;
252    } else {
253      solveMode_ = 10 * (solveMode_ / 10);
254      // redo matrix
255      // last start
256      pivotRow_[numberRows_] = lastStart;
257      preProcess();
258    }
259  }
260  status_ = 0;
261  for (int j = 0; j < numberRows_; j++) {
262    pivotRow_[j + numberDense_] = j;
263  }
264  CoinFactorizationDouble *elements = area;
265  numberGoodU_ = 0;
266  for (int i = 0; i < numberRows_; i++) {
267    int iRow = -1;
268    // Find largest
269    double largest = zeroTolerance_;
270    for (int j = i; j < numberRows_; j++) {
271      double value = fabs(elements[j]);
272      if (value > largest) {
273        largest = value;
274        iRow = j;
275      }
276    }
277    if (iRow >= 0) {
278      if (iRow != i) {
279        // swap
280        assert(iRow > i);
281        CoinFactorizationDouble *elementsA = area;
282        for (int k = 0; k <= i; k++) {
283          // swap
284          CoinFactorizationDouble value = elementsA[i];
285          elementsA[i] = elementsA[iRow];
286          elementsA[iRow] = value;
287          elementsA += numberDense_;
288        }
289        int iPivot = pivotRow_[i + numberDense_];
290        pivotRow_[i + numberDense_] = pivotRow_[iRow + numberDense_];
291        pivotRow_[iRow + numberDense_] = iPivot;
292      }
293      CoinFactorizationDouble pivotValue = 1.0 / elements[i];
294      elements[i] = pivotValue;
295      for (int j = i + 1; j < numberRows_; j++) {
296        elements[j] *= pivotValue;
297      }
298      // Update rest
299      CoinFactorizationDouble *elementsA = elements;
300      for (int k = i + 1; k < numberRows_; k++) {
301        elementsA += numberDense_;
302        // swap
303        if (iRow != i) {
304          CoinFactorizationDouble value = elementsA[i];
305          elementsA[i] = elementsA[iRow];
306          elementsA[iRow] = value;
307        }
308        CoinFactorizationDouble value = elementsA[i];
309        for (int j = i + 1; j < numberRows_; j++) {
310          elementsA[j] -= value * elements[j];
311        }
312      }
313    } else {
314      status_ = -1;
315      break;
316    }
317    numberGoodU_++;
318    elements += numberDense_;
319  }
320  for (int j = 0; j < numberRows_; j++) {
321    int k = pivotRow_[j + numberDense_];
322    pivotRow_[k] = j;
323  }
324  //assert (status_);
325  //if (!status_)
326  //solveMode_=10; // for next time
327  //for (int j=0;j<numberRows_;j++) {
328  //int k = pivotRow_[j+numberDense_];
329  //pivotRow_[k]=j;
330  //}
331  return status_;
332}
333// Makes a non-singular basis by replacing variables
334void CoinAbcDenseFactorization::makeNonSingular(int *sequence)
335{
336  // Replace bad ones by correct slack
337  int *workArea = reinterpret_cast< int * >(workArea_);
338  int i;
339  for (i = 0; i < numberRows_; i++)
340    workArea[i] = -1;
341  for (i = 0; i < numberGoodU_; i++) {
342    int iOriginal = pivotRow_[i + numberDense_];
343    workArea[iOriginal] = i;
344    //workArea[i]=iOriginal;
345  }
346  int lastRow = -1;
347  for (i = 0; i < numberRows_; i++) {
348    if (workArea[i] == -1) {
349      lastRow = i;
350      break;
351    }
352  }
353  assert(lastRow >= 0);
354  for (i = numberGoodU_; i < numberRows_; i++) {
355    assert(lastRow < numberRows_);
356    // Put slack in basis
357    sequence[i] = lastRow;
358    lastRow++;
359    for (; lastRow < numberRows_; lastRow++) {
360      if (workArea[lastRow] == -1)
361        break;
362    }
363  }
364}
365//#define DENSE_PERMUTE
366// Does post processing on valid factorization - putting variables on correct rows
367void CoinAbcDenseFactorization::postProcess(const int *sequence, int *pivotVariable)
368{
369  if ((solveMode_ % 10) == 0) {
370    for (int i = 0; i < numberRows_; i++) {
371      int k = sequence[i];
372#ifdef DENSE_PERMUTE
373      pivotVariable[pivotRow_[i + numberDense_]] = k;
374#else
375      //pivotVariable[pivotRow_[i]]=k;
376      //pivotVariable[pivotRow_[i]]=k;
377      pivotVariable[i] = k;
378      k = pivotRow_[i];
379      pivotRow_[i] = pivotRow_[i + numberDense_];
380      pivotRow_[i + numberDense_] = k;
381#endif
382    }
383  } else {
384    // lapack
385    for (int i = 0; i < numberRows_; i++) {
386      int k = sequence[i];
387      pivotVariable[i] = k;
388    }
389  }
390}
391/* Replaces one Column to basis,
392   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
393   partial update already in U */
394int CoinAbcDenseFactorization::replaceColumn(CoinIndexedVector *regionSparse,
395  int pivotRow,
396  double pivotCheck,
397  bool /*skipBtranU*/,
398  double /*acceptablePivot*/)
399{
400  if (numberPivots_ == maximumPivots_)
401    return 3;
402  CoinFactorizationDouble *elements = elements_ + numberDense_ * numberPivots_;
403  double *region = regionSparse->denseVector();
404  int *regionIndex = regionSparse->getIndices();
405  int numberNonZero = regionSparse->getNumElements();
406  int i;
407  memset(elements, 0, numberRows_ * sizeof(CoinFactorizationDouble));
408  CoinFactorizationDouble pivotValue = pivotCheck;
409  if (fabs(pivotValue) < zeroTolerance_)
410    return 2;
411  pivotValue = 1.0 / pivotValue;
412  if ((solveMode_ % 10) == 0) {
413    if (regionSparse->packedMode()) {
414      for (i = 0; i < numberNonZero; i++) {
415        int iRow = regionIndex[i];
416        double value = region[i];
417#ifdef DENSE_PERMUTE
418        iRow = pivotRow_[iRow]; // permute
419#endif
420        elements[iRow] = value;
421        ;
422      }
423    } else {
424      // not packed! - from user pivot?
425      for (i = 0; i < numberNonZero; i++) {
426        int iRow = regionIndex[i];
427        double value = region[iRow];
428#ifdef DENSE_PERMUTE
429        iRow = pivotRow_[iRow]; // permute
430#endif
431        elements[iRow] = value;
432        ;
433      }
434    }
435    int realPivotRow = pivotRow_[pivotRow];
436    elements[realPivotRow] = pivotValue;
437    pivotRow_[2 * numberDense_ + numberPivots_] = realPivotRow;
438  } else {
439    // lapack
440    if (regionSparse->packedMode()) {
441      for (i = 0; i < numberNonZero; i++) {
442        int iRow = regionIndex[i];
443        double value = region[i];
444        elements[iRow] = value;
445        ;
446      }
447    } else {
448      // not packed! - from user pivot?
449      for (i = 0; i < numberNonZero; i++) {
450        int iRow = regionIndex[i];
451        double value = region[iRow];
452        elements[iRow] = value;
453        ;
454      }
455    }
456    elements[pivotRow] = pivotValue;
457    pivotRow_[2 * numberDense_ + numberPivots_] = pivotRow;
458  }
459  numberPivots_++;
460  return 0;
461}
462/* Checks if can replace one Column to basis,
463   returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots */
464int CoinAbcDenseFactorization::checkReplacePart2(int pivotRow,
465  double btranAlpha,
466  double ftranAlpha,
467#ifdef ABC_LONG_FACTORIZATION
468  long
469#endif
470  double ftAlpha,
471  double acceptablePivot)
472{
473  if (numberPivots_ == maximumPivots_)
474    return 3;
475  if (fabs(ftranAlpha) < zeroTolerance_)
476    return 2;
477  return 0;
478}
479/* Replaces one Column to basis,
480   partial update already in U */
481void CoinAbcDenseFactorization::replaceColumnPart3(const AbcSimplex *model,
482  CoinIndexedVector *regionSparse,
483  CoinIndexedVector *tableauColumn,
484  int pivotRow,
485#ifdef ABC_LONG_FACTORIZATION
486  long
487#endif
488  double alpha)
489{
490  CoinFactorizationDouble *elements = elements_ + numberDense_ * numberPivots_;
491  double *region = tableauColumn->denseVector();
492  int *regionIndex = tableauColumn->getIndices();
493  int numberNonZero = tableauColumn->getNumElements();
494  int i;
495  memset(elements, 0, numberRows_ * sizeof(CoinFactorizationDouble));
496  double pivotValue = 1.0 / alpha;
497  if ((solveMode_ % 10) == 0) {
498    if (tableauColumn->packedMode()) {
499      for (i = 0; i < numberNonZero; i++) {
500        int iRow = regionIndex[i];
501        double value = region[i];
502#ifdef DENSE_PERMUTE
503        iRow = pivotRow_[iRow]; // permute
504#endif
505        elements[iRow] = value;
506        ;
507      }
508    } else {
509      // not packed! - from user pivot?
510      for (i = 0; i < numberNonZero; i++) {
511        int iRow = regionIndex[i];
512        double value = region[iRow];
513#ifdef DENSE_PERMUTE
514        iRow = pivotRow_[iRow]; // permute
515#endif
516        elements[iRow] = value;
517        ;
518      }
519    }
520    int realPivotRow = pivotRow_[pivotRow];
521    elements[realPivotRow] = pivotValue;
522    pivotRow_[2 * numberDense_ + numberPivots_] = realPivotRow;
523  } else {
524    // lapack
525    if (tableauColumn->packedMode()) {
526      for (i = 0; i < numberNonZero; i++) {
527        int iRow = regionIndex[i];
528        double value = region[i];
529        elements[iRow] = value;
530        ;
531      }
532    } else {
533      // not packed! - from user pivot?
534      for (i = 0; i < numberNonZero; i++) {
535        int iRow = regionIndex[i];
536        double value = region[iRow];
537        elements[iRow] = value;
538        ;
539      }
540    }
541    elements[pivotRow] = pivotValue;
542    pivotRow_[2 * numberDense_ + numberPivots_] = pivotRow;
543  }
544  numberPivots_++;
545}
546static void
547CoinAbcDlaswp1(double *COIN_RESTRICT a, int m, int *ipiv)
548{
549  for (int i = 0; i < m; i++) {
550    int ip = ipiv[i];
551    if (ip != i) {
552      double temp = a[i];
553      a[i] = a[ip];
554      a[ip] = temp;
555    }
556  }
557}
558static void
559CoinAbcDlaswp1Backwards(double *COIN_RESTRICT a, int m, int *ipiv)
560{
561  for (int i = m - 1; i >= 0; i--) {
562    int ip = ipiv[i];
563    if (ip != i) {
564      double temp = a[i];
565      a[i] = a[ip];
566      a[ip] = temp;
567    }
568  }
569}
570/* This version has same effect as above with FTUpdate==false
571   so number returned is always >=0 */
572int CoinAbcDenseFactorization::updateColumn(CoinIndexedVector &regionSparse) const
573{
574  double *region = regionSparse.denseVector();
575  CoinBigIndex put = numberDense_ * numberDense_;
576  CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
577  CoinFactorizationDouble *COIN_RESTRICT elements = area;
578  if ((solveMode_ % 10) == 0) {
579    // base factorization L
580    for (int i = 0; i < numberRows_; i++) {
581      double value = region[i];
582      for (int j = i + 1; j < numberRows_; j++) {
583        region[j] -= value * elements[j];
584      }
585      elements += numberDense_;
586    }
587    elements = area + numberRows_ * numberDense_;
588    // base factorization U
589    for (int i = numberRows_ - 1; i >= 0; i--) {
590      elements -= numberDense_;
591      CoinFactorizationDouble value = region[i] * elements[i];
592      region[i] = value;
593      for (int j = 0; j < i; j++) {
594        region[j] -= value * elements[j];
595      }
596    }
597  } else {
598    CoinAbcDlaswp1(region, numberRows_, pivotRow_ + numberDense_);
599    CoinAbcDgetrs('N', numberDense_, area, region);
600  }
601  // now updates
602  elements = elements_;
603  for (int i = 0; i < numberPivots_; i++) {
604    int iPivot = pivotRow_[i + 2 * numberDense_];
605    CoinFactorizationDouble value = region[iPivot] * elements[iPivot];
606    for (int j = 0; j < numberRows_; j++) {
607      region[j] -= value * elements[j];
608    }
609    region[iPivot] = value;
610    elements += numberDense_;
611  }
612  regionSparse.setNumElements(0);
613  regionSparse.scan(0, numberRows_);
614  return 0;
615}
616
617/* Updates one column for dual steepest edge weights (FTRAN) */
618void CoinAbcDenseFactorization::updateWeights(CoinIndexedVector &regionSparse) const
619{
620  updateColumn(regionSparse);
621}
622
623int CoinAbcDenseFactorization::updateTwoColumnsFT(CoinIndexedVector &regionSparse2,
624  CoinIndexedVector &regionSparse3)
625{
626  updateColumn(regionSparse2);
627  updateColumn(regionSparse3);
628  return 0;
629}
630
631/* Updates one column (BTRAN) from regionSparse2
632   regionSparse starts as zero and is zero at end
633   Note - if regionSparse2 packed on input - will be packed on output
634*/
635int CoinAbcDenseFactorization::updateColumnTranspose(CoinIndexedVector &regionSparse2) const
636{
637  double *region = regionSparse2.denseVector();
638  CoinFactorizationDouble *elements = elements_ + numberDense_ * numberPivots_;
639  // updates
640  for (int i = numberPivots_ - 1; i >= 0; i--) {
641    elements -= numberDense_;
642    int iPivot = pivotRow_[i + 2 * numberDense_];
643    CoinFactorizationDouble value = region[iPivot]; //*elements[iPivot];
644    for (int j = 0; j < iPivot; j++) {
645      value -= region[j] * elements[j];
646    }
647    for (int j = iPivot + 1; j < numberRows_; j++) {
648      value -= region[j] * elements[j];
649    }
650    region[iPivot] = value * elements[iPivot];
651  }
652  CoinBigIndex put = numberDense_ * numberDense_;
653  CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
654  if ((solveMode_ % 10) == 0) {
655    // base factorization U
656    elements = area;
657    for (int i = 0; i < numberRows_; i++) {
658      //CoinFactorizationDouble value = region[i]*elements[i];
659      CoinFactorizationDouble value = region[i];
660      for (int j = 0; j < i; j++) {
661        value -= region[j] * elements[j];
662      }
663      //region[i] = value;
664      region[i] = value * elements[i];
665      elements += numberDense_;
666    }
667    // base factorization L
668    elements = area + numberDense_ * numberRows_;
669    for (int i = numberRows_ - 1; i >= 0; i--) {
670      elements -= numberDense_;
671      CoinFactorizationDouble value = region[i];
672      for (int j = i + 1; j < numberRows_; j++) {
673        value -= region[j] * elements[j];
674      }
675      region[i] = value;
676    }
677  } else {
678    CoinAbcDgetrs('T', numberDense_, area, region);
679    CoinAbcDlaswp1Backwards(region, numberRows_, pivotRow_ + numberDense_);
680  }
681  regionSparse2.setNumElements(0);
682  regionSparse2.scan(0, numberRows_);
683  return 0;
684}
685/* Updates one column (FTRAN) */
686void CoinAbcAnyFactorization::updateColumnCpu(CoinIndexedVector &regionSparse, int whichCpu) const
687{
688  updateColumn(regionSparse);
689}
690/* Updates one column (BTRAN) */
691void CoinAbcAnyFactorization::updateColumnTransposeCpu(CoinIndexedVector &regionSparse, int whichCpu) const
692{
693  updateColumnTranspose(regionSparse);
694}
695// Default constructor
696CoinAbcAnyFactorization::CoinAbcAnyFactorization()
697  : pivotTolerance_(1.0e-1)
698  , minimumPivotTolerance_(1.0e-1)
699  , zeroTolerance_(1.0e-13)
700  , relaxCheck_(1.0)
701  , factorElements_(0)
702  , numberRows_(0)
703  , numberGoodU_(0)
704  , maximumPivots_(200)
705  , numberPivots_(0)
706  , status_(-1)
707  , solveMode_(0)
708{
709}
710// Copy constructor
711CoinAbcAnyFactorization::CoinAbcAnyFactorization(const CoinAbcAnyFactorization &other)
712  : pivotTolerance_(other.pivotTolerance_)
713  , minimumPivotTolerance_(other.minimumPivotTolerance_)
714  , zeroTolerance_(other.zeroTolerance_)
715  , relaxCheck_(other.relaxCheck_)
716  , factorElements_(other.factorElements_)
717  , numberRows_(other.numberRows_)
718  , numberGoodU_(other.numberGoodU_)
719  , maximumPivots_(other.maximumPivots_)
720  , numberPivots_(other.numberPivots_)
721  , status_(other.status_)
722  , solveMode_(other.solveMode_)
723{
724}
725// Destructor
726CoinAbcAnyFactorization::~CoinAbcAnyFactorization()
727{
728}
729// = copy
730CoinAbcAnyFactorization &CoinAbcAnyFactorization::operator=(const CoinAbcAnyFactorization &other)
731{
732  if (this != &other) {
733    pivotTolerance_ = other.pivotTolerance_;
734    minimumPivotTolerance_ = other.minimumPivotTolerance_;
735    zeroTolerance_ = other.zeroTolerance_;
736    relaxCheck_ = other.relaxCheck_;
737    factorElements_ = other.factorElements_;
738    numberRows_ = other.numberRows_;
739    numberGoodU_ = other.numberGoodU_;
740    maximumPivots_ = other.maximumPivots_;
741    numberPivots_ = other.numberPivots_;
742    status_ = other.status_;
743    solveMode_ = other.solveMode_;
744  }
745  return *this;
746}
747void CoinAbcAnyFactorization::pivotTolerance(double value)
748{
749  if (value > 0.0 && value <= 1.0) {
750    pivotTolerance_ = CoinMax(value, minimumPivotTolerance_);
751  }
752}
753void CoinAbcAnyFactorization::minimumPivotTolerance(double value)
754{
755  if (value > 0.0 && value <= 1.0) {
756    minimumPivotTolerance_ = value;
757  }
758}
759void CoinAbcAnyFactorization::zeroTolerance(double value)
760{
761  if (value > 0.0 && value < 1.0) {
762    zeroTolerance_ = value;
763  }
764}
765void CoinAbcAnyFactorization::maximumPivots(int value)
766{
767  if (value > maximumPivots_) {
768    delete[] pivotRow_;
769    int n = maximumRows_;
770    if ((n & (BLOCKING8 - 1)) != 0)
771      n += (BLOCKING8 - (n & (BLOCKING8 - 1)));
772    pivotRow_ = new int[2 * n + value];
773  }
774  maximumPivots_ = value;
775}
776// Number of entries in each row
777int *CoinAbcAnyFactorization::numberInRow() const
778{
779  return reinterpret_cast< int * >(workArea_);
780}
781// Number of entries in each column
782int *CoinAbcAnyFactorization::numberInColumn() const
783{
784  return (reinterpret_cast< int * >(workArea_)) + numberRows_;
785}
786// Returns array to put basis starts in
787CoinBigIndex *
788CoinAbcAnyFactorization::starts() const
789{
790  return reinterpret_cast< CoinBigIndex * >(pivotRow_);
791}
792// Returns array to put basis elements in
793CoinFactorizationDouble *
794CoinAbcAnyFactorization::elements() const
795{
796  return elements_;
797}
798// Returns pivot row
799int *CoinAbcAnyFactorization::pivotRow() const
800{
801  return pivotRow_;
802}
803// Returns work area
804CoinFactorizationDouble *
805CoinAbcAnyFactorization::workArea() const
806{
807  return workArea_;
808}
809// Returns int work area
810int *CoinAbcAnyFactorization::intWorkArea() const
811{
812  return reinterpret_cast< int * >(workArea_);
813}
814// Returns permute back
815int *CoinAbcAnyFactorization::permuteBack() const
816{
817  return pivotRow_ + numberRows_;
818}
819// Returns pivot column
820int *CoinAbcAnyFactorization::pivotColumn() const
821{
822  return permute();
823}
824// Returns true if wants tableauColumn in replaceColumn
825bool CoinAbcAnyFactorization::wantsTableauColumn() const
826{
827  return true;
828}
829/* Useful information for factorization
830   0 - iteration number
831   whereFrom is 0 for factorize and 1 for replaceColumn
832*/
833void CoinAbcAnyFactorization::setUsefulInformation(const int *, int)
834{
835}
836
837/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
838*/
Note: See TracBrowser for help on using the repository browser.