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

Last change on this file since 2030 was 1910, checked in by stefan, 7 years ago
  • add configure option --enable-aboca={1,2,3,4,yes,no}
  • compile Aboca source only if --enable-aboca set (instead of compiling empty source files)
  • fix svn properties
  • Property svn:keywords set to Id
File size: 22.9 KB
Line 
1/* $Id: CoinAbcDenseFactorization.cpp 1910 2013-01-27 02:00:13Z forrest $ */
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  if (this != &other) {   
92    gutsOfDestructor();
93    gutsOfInitialize();
94    gutsOfCopy(other);
95  }
96  return *this;
97}
98#define WORK_MULT 2
99void CoinAbcDenseFactorization::gutsOfCopy(const CoinAbcDenseFactorization &other)
100{
101  pivotTolerance_ = other.pivotTolerance_;
102  minimumPivotTolerance_ = other.minimumPivotTolerance_;
103  zeroTolerance_ = other.zeroTolerance_;
104  areaFactor_=other.areaFactor_;
105  numberDense_=other.numberDense_;
106  relaxCheck_ = other.relaxCheck_;
107  numberRows_ = other.numberRows_;
108  maximumRows_ = other.maximumRows_;
109#if ABC_PARALLEL==2
110  parallelMode_=other.parallelMode_;
111#endif
112  maximumSpace_ = other.maximumSpace_;
113  maximumRowsAdjusted_ = other.maximumRowsAdjusted_;
114  solveMode_ = other.solveMode_;
115  numberGoodU_ = other.numberGoodU_;
116  maximumPivots_ = other.maximumPivots_;
117  numberPivots_ = other.numberPivots_;
118  factorElements_ = other.factorElements_;
119  status_ = other.status_;
120  numberSlacks_ = other.numberSlacks_;
121  if (other.pivotRow_) {
122    pivotRow_ = new int [2*maximumRowsAdjusted_+maximumPivots_];
123    CoinMemcpyN(other.pivotRow_,(2*maximumRowsAdjusted_+numberPivots_),pivotRow_);
124    elements_ = new CoinFactorizationDouble [maximumSpace_];
125    CoinMemcpyN(other.elements_,(maximumRowsAdjusted_+numberPivots_)*maximumRowsAdjusted_,elements_);
126    workArea_ = new CoinFactorizationDouble [maximumRowsAdjusted_*WORK_MULT];
127    CoinZeroN(workArea_,maximumRowsAdjusted_*WORK_MULT);
128  } else {
129    elements_ = NULL;
130    pivotRow_ = NULL;
131    workArea_ = NULL;
132  }
133}
134
135//  getAreas.  Gets space for a factorization
136//called by constructors
137void
138CoinAbcDenseFactorization::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
168CoinAbcDenseFactorization::preProcess ()
169{
170  CoinBigIndex put = numberDense_*numberDense_;
171  CoinFactorizationDouble * COIN_RESTRICT area = elements_+maximumSpace_-put;
172  CoinAbcMemset0(area,put);
173  int *indexRow = reinterpret_cast<int *> (elements_+numberRows_*numberRows_);
174  CoinBigIndex * starts = reinterpret_cast<CoinBigIndex *> (pivotRow_); 
175  CoinFactorizationDouble * COIN_RESTRICT column = area;
176  //if (solveMode_==10)
177  //solveMode_=1;
178  if ((solveMode_%10)!=0) {
179    for (int i=0;i<numberRows_;i++) {
180      for (CoinBigIndex j=starts[i];j<starts[i+1];j++) {
181        int iRow = indexRow[j];
182        int iBlock=iRow>>3;
183        iRow=iRow&(BLOCKING8-1);
184        column[iRow+BLOCKING8X8*iBlock]=elements_[j];
185      }
186      if ((i&(BLOCKING8-1))!=(BLOCKING8-1))
187        column += BLOCKING8;
188      else
189        column += numberDense_*BLOCKING8-(BLOCKING8-1)*BLOCKING8;
190    }
191    for (int i=numberRows_;i<numberDense_;i++) {
192      int iRow = i;
193      int iBlock=iRow>>3;
194      iRow=iRow&(BLOCKING8-1);
195      column[iRow+BLOCKING8X8*iBlock]=1.0;
196      if ((i&(BLOCKING8-1))!=(BLOCKING8-1))
197        column += BLOCKING8;
198      //else
199      //column += numberDense_*BLOCKING8-(BLOCKING8-1)*BLOCKING8;
200    }
201  } else {
202    // first or bad
203    for (int i=0;i<numberRows_;i++) {
204      for (CoinBigIndex j=starts[i];j<starts[i+1];j++) {
205        int iRow = indexRow[j];
206        column[iRow]=elements_[j];
207      }
208      column += numberDense_;
209    }
210    for (int i=numberRows_;i<numberDense_;i++) {
211      column[i]=1.0;
212      column += numberDense_;
213    }
214  }
215}
216
217//Does factorization
218int
219CoinAbcDenseFactorization::factor (AbcSimplex * /*model*/)
220{
221  numberPivots_=0;
222  // ? take out
223  //printf("Debug non lapack dense factor\n");
224  //solveMode_&=~1;
225  CoinBigIndex put = numberDense_*numberDense_;
226  CoinFactorizationDouble * COIN_RESTRICT area = elements_+maximumSpace_-put;
227  if ((solveMode_%10)!=0) {
228    // save last start
229    CoinBigIndex lastStart=pivotRow_[numberRows_];
230    status_=CoinAbcDgetrf(numberDense_,numberDense_,area,numberDense_,pivotRow_+numberDense_
231#if ABC_PARALLEL==2
232                          ,parallelMode_
233#endif
234                          );
235    // need to check size of pivots
236    if(!status_) {
237      // OK
238      solveMode_=1+10*(solveMode_/10);
239      numberGoodU_=numberRows_;
240      CoinZeroN(workArea_,2*numberRows_);
241#if 0 //ndef NDEBUG
242      const CoinFactorizationDouble * column = area;
243      double smallest=COIN_DBL_MAX;
244      for (int i=0;i<numberRows_;i++) {
245        if (fabs(column[i])<smallest)
246          smallest = fabs(column[i]);
247        column += numberRows_;
248      }
249      if (smallest<1.0e-8)
250        printf("small el %g\n",smallest);
251#endif
252      return 0;
253    } else {
254      solveMode_=10*(solveMode_/10);
255      // redo matrix
256      // last start
257      pivotRow_[numberRows_]=lastStart;
258      preProcess();
259    }
260  }
261  status_=0;
262  for (int j=0;j<numberRows_;j++) {
263    pivotRow_[j+numberDense_]=j;
264  }
265  CoinFactorizationDouble * elements = area;
266  numberGoodU_=0;
267  for (int i=0;i<numberRows_;i++) {
268    int iRow = -1;
269    // Find largest
270    double largest=zeroTolerance_;
271    for (int j=i;j<numberRows_;j++) {
272      double value = fabs(elements[j]);
273      if (value>largest) {
274        largest=value;
275        iRow=j;
276      }
277    }
278    if (iRow>=0) {
279      if (iRow!=i) {
280        // swap
281        assert (iRow>i);
282        CoinFactorizationDouble * elementsA = area;
283        for (int k=0;k<=i;k++) {
284          // swap
285          CoinFactorizationDouble value = elementsA[i];
286          elementsA[i]=elementsA[iRow];
287          elementsA[iRow]=value;
288          elementsA += numberDense_;
289        }
290        int iPivot = pivotRow_[i+numberDense_];
291        pivotRow_[i+numberDense_]=pivotRow_[iRow+numberDense_];
292        pivotRow_[iRow+numberDense_]=iPivot;
293      }
294      CoinFactorizationDouble pivotValue = 1.0/elements[i];
295      elements[i]=pivotValue;
296      for (int j=i+1;j<numberRows_;j++) {
297        elements[j] *= pivotValue;
298      }
299      // Update rest
300      CoinFactorizationDouble * elementsA = elements;
301      for (int k=i+1;k<numberRows_;k++) {
302        elementsA += numberDense_;
303        // swap
304        if (iRow!=i) {
305          CoinFactorizationDouble value = elementsA[i];
306          elementsA[i]=elementsA[iRow];
307          elementsA[iRow]=value;
308        }
309        CoinFactorizationDouble value = elementsA[i];
310        for (int j=i+1;j<numberRows_;j++) {
311          elementsA[j] -= value * elements[j];
312        }
313      }
314    } else {
315      status_=-1;
316      break;
317    }
318    numberGoodU_++;
319    elements += numberDense_;
320  }
321  for (int j=0;j<numberRows_;j++) {
322    int k = pivotRow_[j+numberDense_];
323    pivotRow_[k]=j;
324  }
325  //assert (status_);
326  //if (!status_)
327  //solveMode_=10; // for next time
328  //for (int j=0;j<numberRows_;j++) {
329  //int k = pivotRow_[j+numberDense_];
330  //pivotRow_[k]=j;
331  //}
332  return status_;
333}
334// Makes a non-singular basis by replacing variables
335void 
336CoinAbcDenseFactorization::makeNonSingular(int * sequence)
337{
338  // Replace bad ones by correct slack
339  int * workArea = reinterpret_cast<int *> (workArea_);
340  int i;
341  for ( i=0;i<numberRows_;i++) 
342    workArea[i]=-1;
343  for ( i=0;i<numberGoodU_;i++) {
344    int iOriginal = pivotRow_[i+numberDense_];
345    workArea[iOriginal]=i;
346    //workArea[i]=iOriginal;
347  }
348  int lastRow=-1;
349  for ( i=0;i<numberRows_;i++) {
350    if (workArea[i]==-1) {
351      lastRow=i;
352      break;
353    }
354  }
355  assert (lastRow>=0);
356  for ( i=numberGoodU_;i<numberRows_;i++) {
357    assert (lastRow<numberRows_);
358    // Put slack in basis
359    sequence[i]=lastRow;
360    lastRow++;
361    for (;lastRow<numberRows_;lastRow++) {
362      if (workArea[lastRow]==-1)
363        break;
364    }
365  }
366}
367//#define DENSE_PERMUTE
368// Does post processing on valid factorization - putting variables on correct rows
369void 
370CoinAbcDenseFactorization::postProcess(const int * sequence, int * pivotVariable)
371{
372  if ((solveMode_%10)==0) {
373    for (int i=0;i<numberRows_;i++) {
374      int k = sequence[i];
375#ifdef DENSE_PERMUTE
376      pivotVariable[pivotRow_[i+numberDense_]]=k;
377#else
378      //pivotVariable[pivotRow_[i]]=k;
379      //pivotVariable[pivotRow_[i]]=k;
380      pivotVariable[i]=k;
381      k=pivotRow_[i];
382      pivotRow_[i] = pivotRow_[i+numberDense_];
383      pivotRow_[i+numberDense_]=k;
384#endif
385    }
386  } else {
387    // lapack
388    for (int i=0;i<numberRows_;i++) {
389      int k = sequence[i];
390      pivotVariable[i]=k;
391    }
392  }
393}
394/* Replaces one Column to basis,
395   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
396   partial update already in U */
397int 
398CoinAbcDenseFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
399                                        int pivotRow,
400                                        double pivotCheck ,
401                                        bool /*skipBtranU*/,
402                                       double /*acceptablePivot*/)
403{
404  if (numberPivots_==maximumPivots_)
405    return 3;
406  CoinFactorizationDouble * elements = elements_ + numberDense_*numberPivots_;
407  double *region = regionSparse->denseVector (  );
408  int *regionIndex = regionSparse->getIndices (  );
409  int numberNonZero = regionSparse->getNumElements (  );
410  int i;
411  memset(elements,0,numberRows_*sizeof(CoinFactorizationDouble));
412  CoinFactorizationDouble pivotValue = pivotCheck;
413  if (fabs(pivotValue)<zeroTolerance_)
414    return 2;
415  pivotValue = 1.0/pivotValue;
416  if ((solveMode_%10)==0) {
417    if (regionSparse->packedMode()) {
418      for (i=0;i<numberNonZero;i++) {
419        int iRow = regionIndex[i];
420        double value = region[i];
421#ifdef DENSE_PERMUTE
422        iRow = pivotRow_[iRow]; // permute
423#endif
424        elements[iRow] = value;;
425      }
426    } else {
427      // not packed! - from user pivot?
428      for (i=0;i<numberNonZero;i++) {
429        int iRow = regionIndex[i];
430        double value = region[iRow];
431#ifdef DENSE_PERMUTE
432        iRow = pivotRow_[iRow]; // permute
433#endif
434        elements[iRow] = value;;
435      }
436    }
437    int realPivotRow = pivotRow_[pivotRow];
438    elements[realPivotRow]=pivotValue;
439    pivotRow_[2*numberDense_+numberPivots_]=realPivotRow;
440  } else {
441    // lapack
442    if (regionSparse->packedMode()) {
443      for (i=0;i<numberNonZero;i++) {
444        int iRow = regionIndex[i];
445        double value = region[i];
446        elements[iRow] = value;;
447      }
448    } else {
449      // not packed! - from user pivot?
450      for (i=0;i<numberNonZero;i++) {
451        int iRow = regionIndex[i];
452        double value = region[iRow];
453        elements[iRow] = value;;
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 
465CoinAbcDenseFactorization::checkReplacePart2 ( int pivotRow,
466                                               double btranAlpha, 
467                                               double ftranAlpha, 
468#ifdef ABC_LONG_FACTORIZATION
469                                               long
470#endif
471                                               double ftAlpha,
472                                               double acceptablePivot)
473{
474  if (numberPivots_==maximumPivots_)
475    return 3;
476  if (fabs(ftranAlpha)<zeroTolerance_)
477    return 2;
478  return 0;
479}
480/* Replaces one Column to basis,
481   partial update already in U */
482void 
483CoinAbcDenseFactorization::replaceColumnPart3 ( const AbcSimplex * model,
484                                           CoinIndexedVector * regionSparse,
485                                           CoinIndexedVector * tableauColumn,
486                                           int pivotRow,
487#ifdef ABC_LONG_FACTORIZATION
488                                                long
489#endif
490                                           double alpha )
491{
492  CoinFactorizationDouble * elements = elements_ + numberDense_*numberPivots_;
493  double *region = tableauColumn->denseVector (  );
494  int *regionIndex = tableauColumn->getIndices (  );
495  int numberNonZero = tableauColumn->getNumElements (  );
496  int i;
497  memset(elements,0,numberRows_*sizeof(CoinFactorizationDouble));
498  double pivotValue = 1.0/alpha;
499  if ((solveMode_%10)==0) {
500    if (tableauColumn->packedMode()) {
501      for (i=0;i<numberNonZero;i++) {
502        int iRow = regionIndex[i];
503        double value = region[i];
504#ifdef DENSE_PERMUTE
505        iRow = pivotRow_[iRow]; // permute
506#endif
507        elements[iRow] = value;;
508      }
509    } else {
510      // not packed! - from user pivot?
511      for (i=0;i<numberNonZero;i++) {
512        int iRow = regionIndex[i];
513        double value = region[iRow];
514#ifdef DENSE_PERMUTE
515        iRow = pivotRow_[iRow]; // permute
516#endif
517        elements[iRow] = value;;
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    } else {
532      // not packed! - from user pivot?
533      for (i=0;i<numberNonZero;i++) {
534        int iRow = regionIndex[i];
535        double value = region[iRow];
536        elements[iRow] = value;;
537      }
538    }
539    elements[pivotRow]=pivotValue;
540    pivotRow_[2*numberDense_+numberPivots_]=pivotRow;
541  }
542  numberPivots_++;
543}
544static void 
545CoinAbcDlaswp1(double * COIN_RESTRICT a,int m, int * ipiv)
546{
547  for (int i=0;i<m;i++) {
548    int ip=ipiv[i];
549    if (ip!=i) {
550      double temp=a[i];
551      a[i]=a[ip];
552      a[ip]=temp;
553    }
554  }
555}
556static void 
557CoinAbcDlaswp1Backwards(double * COIN_RESTRICT a,int m, int * ipiv)
558{
559  for (int i=m-1;i>=0;i--) {
560    int ip=ipiv[i];
561    if (ip!=i) {
562      double temp=a[i];
563      a[i]=a[ip];
564      a[ip]=temp;
565    }
566  }
567}
568/* This version has same effect as above with FTUpdate==false
569   so number returned is always >=0 */
570int 
571CoinAbcDenseFactorization::updateColumn (CoinIndexedVector & regionSparse) const
572{
573  double *region = regionSparse.denseVector (  );
574  CoinBigIndex put = numberDense_*numberDense_;
575  CoinFactorizationDouble * COIN_RESTRICT area = elements_+maximumSpace_-put;
576  CoinFactorizationDouble * COIN_RESTRICT elements = area;
577  if ((solveMode_%10)==0) {
578    // base factorization L
579    for (int i=0;i<numberRows_;i++) {
580      double value = region[i];
581      for (int j=i+1;j<numberRows_;j++) {
582        region[j] -= value*elements[j];
583      }
584      elements += numberDense_;
585    }
586    elements = area+numberRows_*numberDense_;
587    // base factorization U
588    for (int i=numberRows_-1;i>=0;i--) {
589      elements -= numberDense_;
590      CoinFactorizationDouble value = region[i]*elements[i];
591      region[i] = value;
592      for (int j=0;j<i;j++) {
593        region[j] -= value*elements[j];
594      }
595    }
596  } else {
597    CoinAbcDlaswp1(region,numberRows_,pivotRow_+numberDense_);
598    CoinAbcDgetrs('N',numberDense_,area,region);
599  }
600  // now updates
601  elements = elements_;
602  for (int i=0;i<numberPivots_;i++) {
603    int iPivot = pivotRow_[i+2*numberDense_];
604    CoinFactorizationDouble value = region[iPivot]*elements[iPivot];
605    for (int j=0;j<numberRows_;j++) {
606      region[j] -= value*elements[j];
607    }
608    region[iPivot] = value;
609    elements += numberDense_;
610  }
611  regionSparse.setNumElements(0);
612  regionSparse.scan(0,numberRows_);
613  return 0;
614}
615
616/* Updates one column for dual steepest edge weights (FTRAN) */
617void 
618CoinAbcDenseFactorization::updateWeights ( CoinIndexedVector & regionSparse) const
619{
620  updateColumn(regionSparse);
621}
622
623int 
624CoinAbcDenseFactorization::updateTwoColumnsFT(CoinIndexedVector & regionSparse2,
625                                              CoinIndexedVector & regionSparse3)
626{
627  updateColumn(regionSparse2);
628  updateColumn(regionSparse3);
629  return 0;
630}
631
632/* Updates one column (BTRAN) from regionSparse2
633   regionSparse starts as zero and is zero at end
634   Note - if regionSparse2 packed on input - will be packed on output
635*/
636int 
637CoinAbcDenseFactorization::updateColumnTranspose ( CoinIndexedVector & regionSparse2) const
638{
639  double *region = regionSparse2.denseVector (  );
640  CoinFactorizationDouble * elements = elements_+numberDense_*numberPivots_;
641  // updates
642  for (int i=numberPivots_-1;i>=0;i--) {
643    elements -= numberDense_;
644    int iPivot = pivotRow_[i+2*numberDense_];
645    CoinFactorizationDouble value = region[iPivot]; //*elements[iPivot];
646    for (int j=0;j<iPivot;j++) {
647      value -= region[j]*elements[j];
648    }
649    for (int j=iPivot+1;j<numberRows_;j++) {
650      value -= region[j]*elements[j];
651    }
652    region[iPivot] = value*elements[iPivot];
653  }
654  CoinBigIndex put = numberDense_*numberDense_;
655  CoinFactorizationDouble * COIN_RESTRICT area = elements_+maximumSpace_-put;
656  if ((solveMode_%10)==0) {
657    // base factorization U
658    elements = area;
659    for (int i=0;i<numberRows_;i++) {
660      //CoinFactorizationDouble value = region[i]*elements[i];
661      CoinFactorizationDouble value = region[i];
662      for (int j=0;j<i;j++) {
663        value -= region[j]*elements[j];
664      }
665      //region[i] = value;
666      region[i] = value*elements[i];
667      elements += numberDense_;
668    }
669    // base factorization L
670    elements = area+numberDense_*numberRows_;
671    for (int i=numberRows_-1;i>=0;i--) {
672      elements -= numberDense_;
673      CoinFactorizationDouble value = region[i];
674      for (int j=i+1;j<numberRows_;j++) {
675        value -= region[j]*elements[j];
676      }
677      region[i] = value;
678    }
679  } else {
680    CoinAbcDgetrs('T',numberDense_,area,region);
681    CoinAbcDlaswp1Backwards(region,numberRows_,pivotRow_+numberDense_);
682  }
683  regionSparse2.setNumElements(0);
684  regionSparse2.scan(0,numberRows_);
685  return 0;
686}
687/* Updates one column (FTRAN) */
688void 
689CoinAbcAnyFactorization::updateColumnCpu ( CoinIndexedVector & regionSparse,int whichCpu) const
690{ updateColumn(regionSparse);}
691/* Updates one column (BTRAN) */
692void 
693CoinAbcAnyFactorization::updateColumnTransposeCpu ( CoinIndexedVector & regionSparse,int whichCpu) const
694{ updateColumnTranspose(regionSparse);}
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 
766CoinAbcAnyFactorization::maximumPivots (  int value )
767{
768  if (value>maximumPivots_) {
769    delete [] pivotRow_;
770    int n=maximumRows_;
771    if ((n&(BLOCKING8-1))!=0)
772    n += (BLOCKING8-(n&(BLOCKING8-1)));
773    pivotRow_ = new int[2*n+value];
774  }
775  maximumPivots_ = value;
776}
777// Number of entries in each row
778int * 
779CoinAbcAnyFactorization::numberInRow() const
780{ return reinterpret_cast<int *> (workArea_);}
781// Number of entries in each column
782int * 
783CoinAbcAnyFactorization::numberInColumn() const
784{ return (reinterpret_cast<int *> (workArea_))+numberRows_;}
785// Returns array to put basis starts in
786CoinBigIndex * 
787CoinAbcAnyFactorization::starts() const
788{ return reinterpret_cast<CoinBigIndex *> (pivotRow_);}
789// Returns array to put basis elements in
790CoinFactorizationDouble * 
791CoinAbcAnyFactorization::elements() const
792{ return elements_;}
793// Returns pivot row
794int * 
795CoinAbcAnyFactorization::pivotRow() const
796{ return pivotRow_;}
797// Returns work area
798CoinFactorizationDouble * 
799CoinAbcAnyFactorization::workArea() const
800{ return workArea_;}
801// Returns int work area
802int * 
803CoinAbcAnyFactorization::intWorkArea() const
804{ return reinterpret_cast<int *> (workArea_);}
805// Returns permute back
806int * 
807CoinAbcAnyFactorization::permuteBack() const
808{ return pivotRow_+numberRows_;}
809// Returns pivot column
810int * 
811CoinAbcAnyFactorization::pivotColumn() const
812{ return permute();}
813// Returns true if wants tableauColumn in replaceColumn
814bool
815CoinAbcAnyFactorization::wantsTableauColumn() const
816{ return true;}
817/* Useful information for factorization
818   0 - iteration number
819   whereFrom is 0 for factorize and 1 for replaceColumn
820*/
821void 
822CoinAbcAnyFactorization::setUsefulInformation(const int * ,int )
823{ }
Note: See TracBrowser for help on using the repository browser.