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

Last change on this file since 1879 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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