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

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

minor changes to implement Aboca

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