source: trunk/Clp/src/CoinAbcBaseFactorization5.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: 26.8 KB
Line 
1/* $Id: CoinAbcBaseFactorization5.cpp 1910 2013-01-27 02:00:13Z forrest $ */
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
222CoinAbcTypeFactorization::updateColumnR ( CoinIndexedVector * regionSparse
223#if ABC_SMALL<2
224                                      , CoinAbcStatistics & statistics
225#endif
226#if ABC_PARALLEL
227                                          ,int whichSparse
228#endif
229                                      ) const
230{
231
232  if ( numberR_ ) {
233    CoinSimplexDouble * COIN_RESTRICT region = regionSparse->denseVector (  );
234#if ABC_SMALL<3
235    CoinSimplexInt * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
236    CoinSimplexInt numberNonZero = regionSparse->getNumElements (  );
237#endif
238   
239    const CoinBigIndex *  COIN_RESTRICT startColumn = startColumnRAddress_-numberRows_;
240    const CoinSimplexInt *  COIN_RESTRICT indexRow = indexRowRAddress_;
241    CoinSimplexInt * COIN_RESTRICT indexRowR2 = indexRowRAddress_+lengthAreaR_;
242    if (gotRCopy()) 
243      indexRowR2 += lengthAreaR_;
244    const CoinFactorizationDouble *  COIN_RESTRICT element = elementRAddress_;
245    const CoinSimplexInt *  COIN_RESTRICT permute = permuteAddress_;
246    instrument_do("CoinAbcFactorizationUpdateR",sizeR);
247#if ABC_SMALL<2
248    // Size of R
249    CoinSimplexDouble sizeR=startColumnRAddress_[numberR_];
250    // Work out very dubious idea of what would be fastest
251    // Average
252    CoinSimplexDouble averageR = sizeR/(static_cast<CoinSimplexDouble> (numberRowsExtra_));
253    // weights (relative to actual work)
254    CoinSimplexDouble setMark = 0.1; // setting mark
255    CoinSimplexDouble test1= 1.0; // starting ftran (without testPivot)
256    CoinSimplexDouble testPivot = 2.0; // Seeing if zero etc
257    CoinSimplexDouble startDot=2.0; // For starting dot product version
258    // For final scan
259    CoinSimplexDouble final = numberNonZero*1.0;
260    CoinSimplexDouble methodTime0;
261    // For first type
262    methodTime0 = numberPivots_ * (testPivot + ((static_cast<double> (numberNonZero))/(static_cast<double> (numberRows_))
263                                                * 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      {
286        instrument_start("CoinAbcFactorizationUpdateRSparse1",numberRows_);
287        // mark known to be zero
288        CoinCheckZero * COIN_RESTRICT mark = markRowAddress_;
289#if ABC_PARALLEL
290        if (whichSparse) {
291          mark=reinterpret_cast<CoinCheckZero *>(reinterpret_cast<char *>(mark)+whichSparse*sizeSparseArray_);
292        }
293#endif
294        // mark all rows which will be permuted
295        for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++ ) {
296          CoinSimplexInt iRow = permute[i];
297          mark[iRow]=1;
298        }
299        // we have another copy of R in R
300        const CoinFactorizationDouble *  COIN_RESTRICT elementR = elementRAddress_ + lengthAreaR_;
301        const CoinSimplexInt *  COIN_RESTRICT indexRowR = indexRowRAddress_ + lengthAreaR_;
302        const CoinBigIndex *  COIN_RESTRICT startR = startColumnRAddress_+maximumPivots_+1;
303        // For current list order does not matter as
304        // only affects end
305        CoinSimplexInt newNumber=0;
306        for (CoinSimplexInt i = 0; i < numberNonZero; i++ ) {
307          CoinSimplexInt iRow = regionIndex[i];
308          CoinFactorizationDouble pivotValue = region[iRow];
309          assert (region[iRow]);
310          if (!mark[iRow]) {
311            regionIndex[newNumber++]=iRow;
312          }
313          CoinSimplexInt kRow=permute[iRow];
314          CoinSimplexInt number = numberInColumnPlus[kRow];
315          instrument_add(number);
316          if (TEST_INT_NONZERO(number)) {
317            CoinBigIndex start=startR[kRow];
318#ifndef INLINE_IT
319            CoinBigIndex end = start+number;
320            for (CoinBigIndex j = start; j < end; j ++ ) {
321              CoinFactorizationDouble value = elementR[j];
322              CoinSimplexInt jRow = indexRowR[j];
323              region[jRow] -= pivotValue*value;
324            }
325#else
326            CoinAbcScatterUpdate(number,pivotValue,elementR+start,indexRowR+start,region);
327#endif
328          }
329        }
330        instrument_start("CoinAbcFactorizationUpdateRSparse2",numberRows_);
331        numberNonZero = newNumber;
332        for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++ ) {
333          //move using permute_ (stored in inverse fashion)
334          CoinSimplexInt iRow = permute[i];
335          CoinFactorizationDouble pivotValue = region[iRow]+region[i];
336          //zero out pre-permuted
337          region[iRow] = 0.0;
338          region[i]=0.0;
339          if ( !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue )  ) {
340            int jRow=pivotRowBack[i];
341            assert (!region[jRow]);
342            region[jRow]=pivotValue;
343            CoinSimplexInt number = numberInColumnPlus[i];
344            instrument_add(number);
345            if (TEST_INT_NONZERO(number)) {
346              CoinBigIndex start=startR[i];
347#ifndef INLINE_IT
348              CoinBigIndex end = start+number;
349              for (CoinBigIndex j = start; j < end; j ++ ) {
350                CoinFactorizationDouble value = elementR[j];
351                CoinSimplexInt jRow = indexRowR[j];
352                region[jRow] -= pivotValue*value;
353              }
354#else
355              CoinAbcScatterUpdate(number,pivotValue,elementR+start,indexRowR+start,region);
356#endif
357            }
358          }
359        }
360        instrument_end();
361        for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++ ) {
362          CoinSimplexInt iRow = permute[i];
363          assert (iRow<numberRows_);
364          if (mark[iRow]) {
365            mark[iRow]=0;
366            CoinExponent expValue=ABC_EXPONENT(region[iRow]);
367            if (expValue) {
368              if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
369                regionIndex[numberNonZero++]=iRow;
370              } else {
371                region[iRow]=0.0;
372              }
373            }
374          }
375        }
376      }
377      break;
378    case 2:
379#endif
380      {
381        instrument_start("CoinAbcFactorizationUpdateRDense",numberRows_);
382        CoinBigIndex start = startColumn[numberRows_];
383        for (CoinSimplexInt i = numberRows_; i < numberRowsExtra_; i++ ) {
384          //move using permute_ (stored in inverse fashion)
385          CoinBigIndex end = startColumn[i+1];
386          CoinSimplexInt iRow = pivotRowBack[i];
387          assert (iRow<numberRows_);
388          CoinFactorizationDouble pivotValue = region[iRow];
389          instrument_add(end-start);
390          if (TEST_INT_NONZERO(end-start)) {
391#ifndef INLINE_IT2
392            for (CoinBigIndex j = start; j < end; j ++ ) {
393              CoinFactorizationDouble value = element[j];
394              CoinSimplexInt jRow = indexRow[j];
395              value *= region[jRow];
396              pivotValue -= value;
397            }
398#else
399            pivotValue+=CoinAbcGatherUpdate(end-start,element+start,indexRow+start,region);
400#endif
401            start=end;
402          }
403#if ABC_SMALL<3
404          if ( !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue ) ) {
405            if (!region[iRow])
406              regionIndex[numberNonZero++] = iRow;
407            region[iRow] = pivotValue;
408          } else {
409            if (region[iRow])
410              region[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
411          }
412#else
413            region[iRow] = pivotValue;
414#endif
415        }
416        instrument_end();
417#if ABC_SMALL<3
418        // pack down
419        CoinSimplexInt n=numberNonZero;
420        numberNonZero=0;
421        for (CoinSimplexInt i=0;i<n;i++) {
422          CoinSimplexInt indexValue = regionIndex[i];
423          CoinExponent expValue=ABC_EXPONENT(region[indexValue]);
424          if (expValue) {
425            if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
426              regionIndex[numberNonZero++]=indexValue; // needed? is region used
427            } else {
428              region[indexValue]=0.0;
429            }
430          }
431        }
432#endif
433      }
434#if ABC_SMALL<2
435      break;
436    }
437#endif
438#if ABC_SMALL<3
439    //set counts
440    regionSparse->setNumElements ( numberNonZero );
441#endif
442  }
443#if ABC_SMALL<2
444  if (factorizationStatistics()) 
445    statistics.countAfterR_ += regionSparse->getNumElements();
446#endif
447}
448#ifdef EARLY_FACTORIZE
449#include "AbcSimplex.hpp"
450#include "AbcMatrix.hpp"
451// Returns -2 if can't, -1 if singular, -99 memory, 0 OK
452int 
453CoinAbcTypeFactorization::factorize (AbcSimplex * model, CoinIndexedVector & stuff)
454{
455  AbcMatrix * matrix = model->abcMatrix();
456  setStatus(-99);
457  int * COIN_RESTRICT pivotTemp = stuff.getIndices();
458  int numberColumnBasic = stuff.getNumElements();
459  // compute how much in basis
460  int numberSlacks = numberRows_ - numberColumnBasic;
461 
462  CoinBigIndex numberElements = numberSlacks+matrix->countBasis(pivotTemp + numberSlacks,
463                                         numberColumnBasic);
464  // Not needed for dense
465  numberElements = 3 * numberRows_ + 3 * numberElements + 20000;
466  getAreas ( numberRows_, numberSlacks + numberColumnBasic, numberElements,
467             2 * numberElements );
468  numberSlacks_=numberSlacks;
469  // Fill in counts so we can skip part of preProcess
470  // This is NOT needed for dense but would be needed for later versions
471  CoinFactorizationDouble *  COIN_RESTRICT elementU;
472  int *  COIN_RESTRICT indexRowU;
473  CoinBigIndex *  COIN_RESTRICT startColumnU;
474  int *  COIN_RESTRICT numberInRow;
475  int *  COIN_RESTRICT numberInColumn;
476  elementU = this->elements();
477  indexRowU = this->indices();
478  startColumnU = this->starts();
479  numberInRow = this->numberInRow();
480  numberInColumn = this->numberInColumn();
481  CoinZeroN ( numberInRow, numberRows_+1  );
482  CoinZeroN ( numberInColumn, numberRows_ );
483  for (int i = 0; i < numberSlacks_; i++) {
484    int iRow = pivotTemp[i];
485    indexRowU[i] = iRow;
486    startColumnU[i] = i;
487    elementU[i] = 1.0;
488    numberInRow[iRow] = 1;
489    numberInColumn[i] = 1;
490  }
491  startColumnU[numberSlacks_] = numberSlacks_;
492  matrix->fillBasis(pivotTemp + numberSlacks_,
493                    numberColumnBasic,
494                    indexRowU,
495                    startColumnU + numberSlacks_,
496                    numberInRow,
497                    numberInColumn + numberSlacks_,
498                    elementU);
499  numberElements = startColumnU[numberRows_-1]
500    + numberInColumn[numberRows_-1];
501  preProcess ( );
502  factor (model);
503  if (status() == 0 ) {
504    // Put sequence numbers in workArea
505    int * savePivots = reinterpret_cast<int *>(workAreaAddress_);
506    CoinAbcMemcpy(savePivots,pivotTemp,numberRows_);
507    postProcess(pivotTemp, savePivots);
508    return 0;
509  } else {
510    return -2;
511  }
512}
513// 0 success, -1 can't +1 accuracy problems
514int 
515CoinAbcTypeFactorization::replaceColumns ( const AbcSimplex * model,
516                                           CoinIndexedVector & stuff,
517                                           int firstPivot,int lastPivot,
518                                           bool cleanUp)
519{
520  // Could skip some if goes in and out (only if easy as then no alpha check)
521  // Put sequence numbers in workArea
522  int * savePivots = reinterpret_cast<int *>(workAreaAddress_);
523  const int * pivotIndices = stuff.getIndices();
524  CoinSimplexDouble * pivotValues = stuff.denseVector();
525  int savePivot = stuff.capacity();
526  savePivot--;
527  savePivot -= 2*firstPivot;
528  int numberDo=lastPivot-firstPivot;
529  // Say clear
530  stuff.setNumElements(0);
531  bool badUpdates=false;
532  for (int iPivot=0;iPivot<numberDo;iPivot++) {
533    int sequenceOut = pivotIndices[--savePivot];
534    CoinSimplexDouble alpha=pivotValues[savePivot];
535    int sequenceIn = pivotIndices[--savePivot];
536    CoinSimplexDouble btranAlpha=pivotValues[savePivot];
537    int pivotRow;
538    for (pivotRow=0;pivotRow<numberRows_;pivotRow++) {
539      if (sequenceOut==savePivots[pivotRow])
540        break;
541    }
542    assert (pivotRow<numberRows_);
543    savePivots[pivotRow]=sequenceIn;
544    model->unpack(stuff,sequenceIn);
545    updateColumnFTPart1(stuff);
546    stuff.clear();
547    //checkReplacePart1a(&stuff,pivotRow);
548    CoinSimplexDouble ftAlpha=checkReplacePart1(&stuff,pivotRow);
549    // may need better check
550    if (checkReplacePart2(pivotRow,btranAlpha,alpha,ftAlpha)>1) {
551      badUpdates=true;
552      printf("Inaccuracy ? btranAlpha %g ftranAlpha %g ftAlpha %g\n",
553             btranAlpha,alpha,ftAlpha);
554      break;
555    }
556    replaceColumnPart3(model,&stuff,NULL,pivotRow,ftAlpha);
557  }
558  int flag;
559  if (!badUpdates) {
560    flag=1;
561    if (cleanUp) {
562      CoinAbcMemcpy(model->pivotVariable(),savePivots,numberRows_);
563    }
564  } else {
565    flag=-1;
566  }
567  stuff.setNumElements(flag);
568  return flag >0 ? 0 : flag;
569}
570#endif
571#ifdef ABC_ORDERED_FACTORIZATION
572// Permute in for Ftran
573void 
574CoinAbcTypeFactorization::permuteInForFtran(CoinIndexedVector & regionSparse,
575                                            bool full) const
576{
577  int numberNonZero=regionSparse.getNumElements();
578  const int * COIN_RESTRICT permuteIn = permuteAddress_+maximumRowsExtra_+1;
579  double * COIN_RESTRICT region = regionSparse.denseVector();
580  double * COIN_RESTRICT tempRegion = region+numberRows_;
581  int * COIN_RESTRICT index = regionSparse.getIndices();
582  if ((numberNonZero<<1)>numberRows_||full) {
583    CoinAbcMemcpy(tempRegion,region,numberRows_);
584    CoinAbcMemset0(region,numberRows_);
585    numberNonZero=0;
586    for (int i=0;i<numberRows_;i++) {
587      double value=tempRegion[i];
588      if (value) {
589        tempRegion[i]=0.0;
590        int iRow=permuteIn[i];
591        region[iRow]=value;
592        index[numberNonZero++]=iRow;
593      }
594    }
595    regionSparse.setNumElements(numberNonZero);
596  } else {
597    for (int i=0;i<numberNonZero;i++) {
598      int iRow=index[i];
599      double value=region[iRow];
600      region[iRow]=0.0;
601      iRow=permuteIn[iRow];
602      tempRegion[iRow]=value;
603      index[i]=iRow;
604    }
605    // and back
606    for (int i=0;i<numberNonZero;i++) {
607      int iRow=index[i];
608      double value=tempRegion[iRow];
609      tempRegion[iRow]=0.0;
610      region[iRow]=value;
611    }
612  }
613}
614// Permute in for Btran and multiply by pivot Region
615void 
616CoinAbcTypeFactorization::permuteInForBtranAndMultiply(CoinIndexedVector & regionSparse,
617                                                       bool full) const
618{
619  int numberNonZero=regionSparse.getNumElements();
620  double * COIN_RESTRICT region = regionSparse.denseVector();
621  double * COIN_RESTRICT tempRegion = region+numberRows_;
622  int * COIN_RESTRICT index = regionSparse.getIndices();
623  const CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegionAddress_;
624  if (full) {
625    numberNonZero=0;
626    for (int iRow=0;iRow<numberRows_;iRow++) {
627      double value = region[iRow];
628      if (value) {
629        region[iRow]=0.0;
630        tempRegion[iRow] = value*pivotRegion[iRow];
631        index[numberNonZero++]=iRow;
632      }
633    }
634    regionSparse.setNumElements(numberNonZero);
635  } else {
636    for (int i=0;i<numberNonZero;i++) {
637      int iRow=index[i];
638      double value = region[iRow];
639      region[iRow]=0.0;
640      tempRegion[iRow] = value*pivotRegion[iRow];
641    }
642  }
643  regionSparse.setDenseVector(tempRegion);
644#if 0
645  if (numberNonZero)
646    printf("BtranIn vector %p region %p regionTemp %p\n",
647         &regionSparse,region,tempRegion);
648  else
649    printf("BtranIn no els! vector %p region %p regionTemp %p\n",
650         &regionSparse,region,tempRegion);
651#endif
652}
653// Permute out for Btran
654void 
655CoinAbcTypeFactorization::permuteOutForBtran(CoinIndexedVector & regionSparse) const
656{
657  int numberNonZero=regionSparse.getNumElements();
658  const int * COIN_RESTRICT permuteIn = permuteAddress_+maximumRowsExtra_+1;
659  const int * COIN_RESTRICT permuteOut = permuteIn+numberRows_;
660  CoinFactorizationDouble * COIN_RESTRICT tempRegion = denseVector(regionSparse);
661  CoinFactorizationDouble * COIN_RESTRICT region = tempRegion-numberRows_;
662  CoinSimplexInt * COIN_RESTRICT index = regionSparse.getIndices();
663  regionSparse.setDenseVector(region);
664  //printf("BtranOut vector %p region %p regionTemp %p\n",
665  //     &regionSparse,region,tempRegion);
666  for (int i=0;i<numberNonZero;i++) {
667    int iRow=index[i];
668    double value=tempRegion[iRow];
669    tempRegion[iRow]=0.0;
670    iRow=permuteOut[iRow];
671    region[iRow]=value;
672    index[i]=iRow;
673  }
674}
675#endif
676#if COIN_BIG_DOUBLE==1
677// To a work array and associate vector
678void 
679CoinAbcTypeFactorization::toLongArray(CoinIndexedVector * vector,int which) const
680{
681  assert (which>=0&&which<FACTOR_CPU);
682  assert (!associatedVector_[which]);
683  associatedVector_[which]=vector;
684  CoinSimplexDouble * COIN_RESTRICT region = vector->denseVector (  );
685  CoinSimplexInt * COIN_RESTRICT regionIndex = vector->getIndices (  );
686  CoinSimplexInt numberNonZero = vector->getNumElements (  );
687  long double * COIN_RESTRICT longRegion = longArray_[which].array();
688  assert (!vector->packedMode());
689  // could check all of longRegion for zero but this should trap most errors
690  for (int i=0;i<numberNonZero;i++) {
691    int iRow=regionIndex[i];
692    assert (!longRegion[iRow]);
693    longRegion[iRow]=region[iRow];
694    region[iRow]=0.0;
695  }
696}
697// From a work array and dis-associate vector
698void 
699CoinAbcTypeFactorization::fromLongArray(CoinIndexedVector * vector) const
700{
701  int which;
702  for (which=0;which<FACTOR_CPU;which++) 
703    if (vector==associatedVector_[which])
704      break;
705  assert (which<FACTOR_CPU);
706  associatedVector_[which]=NULL;
707  CoinSimplexDouble * COIN_RESTRICT region = vector->denseVector (  );
708  CoinSimplexInt * COIN_RESTRICT regionIndex = vector->getIndices (  );
709  CoinSimplexInt numberNonZero = vector->getNumElements (  );
710  long double * COIN_RESTRICT longRegion = longArray_[which].array();
711  assert (!vector->packedMode());
712  // could check all of region for zero but this should trap most errors
713  for (int i=0;i<numberNonZero;i++) {
714    int iRow=regionIndex[i];
715    assert (!region[iRow]);
716    region[iRow]=longRegion[iRow];
717    longRegion[iRow]=0.0;
718  }
719}
720// From a work array and dis-associate vector
721void 
722CoinAbcTypeFactorization::fromLongArray(int which) const
723{
724  assert (which<FACTOR_CPU);
725  CoinIndexedVector * vector = associatedVector_[which];
726  associatedVector_[which]=NULL;
727  CoinSimplexDouble * COIN_RESTRICT region = vector->denseVector (  );
728  CoinSimplexInt * COIN_RESTRICT regionIndex = vector->getIndices (  );
729  CoinSimplexInt numberNonZero = vector->getNumElements (  );
730  long double * COIN_RESTRICT longRegion = longArray_[which].array();
731  assert (!vector->packedMode());
732  // could check all of region for zero but this should trap most errors
733  for (int i=0;i<numberNonZero;i++) {
734    int iRow=regionIndex[i];
735    assert (!region[iRow]);
736    region[iRow]=longRegion[iRow];
737    longRegion[iRow]=0.0;
738  }
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
752long double * 
753CoinAbcTypeFactorization::denseVector(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}
762// Returns long double * associated with vector
763const long double * 
764CoinAbcTypeFactorization::denseVector(const CoinIndexedVector * vector) const
765{
766  int which;
767  for (which=0;which<FACTOR_CPU;which++) 
768    if (vector==associatedVector_[which])
769      break;
770  assert (which<FACTOR_CPU);
771  return longArray_[which].array();
772}
773void
774CoinAbcTypeFactorization::scan(CoinIndexedVector * vector) const
775{
776  int numberNonZero=0;
777  int * COIN_RESTRICT index = vector->getIndices();
778  CoinFactorizationDouble * COIN_RESTRICT region = denseVector(vector);
779  for (int i=0;i<numberRows_;i++) {
780    if (fabs(region[i])>zeroTolerance_)
781      index[numberNonZero++]=i;
782    else
783      region[i]=0.0;
784  }
785  vector->setNumElements(numberNonZero);
786}
787// Clear all hidden arrays
788void 
789CoinAbcTypeFactorization::clearHiddenArrays()
790{
791  for (int which=0;which<FACTOR_CPU;which++) {
792    if (associatedVector_[which]) {
793      CoinIndexedVector * vector = associatedVector_[which];
794      associatedVector_[which]=NULL;
795      CoinFactorizationDouble * COIN_RESTRICT region = longArray_[which].array();
796      CoinSimplexInt * COIN_RESTRICT regionIndex = vector->getIndices (  );
797      CoinSimplexInt numberNonZero = vector->getNumElements (  );
798    }
799  }
800}
801#endif
802#endif
Note: See TracBrowser for help on using the repository browser.