1 /* $Id: CoinAbcBaseFactorization3.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
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 CILK_CONFLICT > 0
17 // for conflicts
18 extern int cilk_conflict;
19 #endif
20 //:class CoinAbcTypeFactorization. Deals with Factorization and Updates
21
22 /* Updates one column for dual steepest edge weights (FTRAN) */
updateWeights(CoinIndexedVector & regionSparse) const23 void CoinAbcTypeFactorization::updateWeights(CoinIndexedVector ®ionSparse) const
24 {
25 toLongArray(®ionSparse, 1);
26 #ifdef ABC_ORDERED_FACTORIZATION
27 // Permute in for Ftran
28 permuteInForFtran(regionSparse);
29 #endif
30 // ******* L
31 updateColumnL(®ionSparse
32 #if ABC_SMALL < 2
33 ,
34 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
35 #endif
36 #if ABC_PARALLEL
37 // can re-use 0 which would have been used for btran
38 ,
39 0
40 #endif
41 );
42 //row bits here
43 updateColumnR(®ionSparse
44 #if ABC_SMALL < 2
45 ,
46 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
47 #endif
48 #if ABC_PARALLEL
49 ,
50 0
51 #endif
52 );
53 //update counts
54 // ******* U
55 updateColumnU(®ionSparse
56 #if ABC_SMALL < 2
57 ,
58 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
59 #endif
60 #if ABC_PARALLEL
61 ,
62 0
63 #endif
64 );
65 #if ABC_DEBUG
66 regionSparse.checkClean();
67 #endif
68 }
updateColumn(CoinIndexedVector & regionSparse) const69 CoinSimplexInt CoinAbcTypeFactorization::updateColumn(CoinIndexedVector ®ionSparse)
70 const
71 {
72 toLongArray(®ionSparse, 0);
73 #ifdef ABC_ORDERED_FACTORIZATION
74 // Permute in for Ftran
75 permuteInForFtran(regionSparse);
76 #endif
77 // ******* L
78 updateColumnL(®ionSparse
79 #if ABC_SMALL < 2
80 ,
81 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
82 #endif
83 );
84 //row bits here
85 updateColumnR(®ionSparse
86 #if ABC_SMALL < 2
87 ,
88 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
89 #endif
90 );
91 //update counts
92 // ******* U
93 updateColumnU(®ionSparse
94 #if ABC_SMALL < 2
95 ,
96 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
97 #endif
98 );
99 #if ABC_DEBUG
100 regionSparse.checkClean();
101 #endif
102 return regionSparse.getNumElements();
103 }
104 /* Updates one full column (FTRAN) */
updateFullColumn(CoinIndexedVector & regionSparse) const105 void CoinAbcTypeFactorization::updateFullColumn(CoinIndexedVector ®ionSparse) const
106 {
107 #ifndef ABC_ORDERED_FACTORIZATION
108 regionSparse.setNumElements(0);
109 regionSparse.scan(0, numberRows_);
110 #else
111 permuteInForFtran(regionSparse, true);
112 #endif
113 if (regionSparse.getNumElements()) {
114 toLongArray(®ionSparse, 1);
115 // ******* L
116 updateColumnL(®ionSparse
117 #if ABC_SMALL < 2
118 ,
119 reinterpret_cast< CoinAbcStatistics & >(ftranFullCountInput_)
120 #endif
121 #if ABC_PARALLEL
122 ,
123 1
124 #endif
125 );
126 //row bits here
127 updateColumnR(®ionSparse
128 #if ABC_SMALL < 2
129 ,
130 reinterpret_cast< CoinAbcStatistics & >(ftranFullCountInput_)
131 #endif
132 #if ABC_PARALLEL
133 ,
134 1
135 #endif
136 );
137 //update counts
138 // ******* U
139 updateColumnU(®ionSparse
140 #if ABC_SMALL < 2
141 ,
142 reinterpret_cast< CoinAbcStatistics & >(ftranFullCountInput_)
143 #endif
144 #if ABC_PARALLEL
145 ,
146 1
147 #endif
148 );
149 fromLongArray(1);
150 #if ABC_DEBUG
151 regionSparse.checkClean();
152 #endif
153 }
154 }
155 // move stuff like this into CoinAbcHelperFunctions.hpp
multiplyIndexed(CoinSimplexInt number,const CoinFactorizationDouble * COIN_RESTRICT multiplier,const CoinSimplexInt * COIN_RESTRICT thisIndex,CoinFactorizationDouble * COIN_RESTRICT region)156 inline void multiplyIndexed(CoinSimplexInt number,
157 const CoinFactorizationDouble *COIN_RESTRICT multiplier,
158 const CoinSimplexInt *COIN_RESTRICT thisIndex,
159 CoinFactorizationDouble *COIN_RESTRICT region)
160 {
161 for (CoinSimplexInt i = 0; i < number; i++) {
162 CoinSimplexInt iRow = thisIndex[i];
163 CoinSimplexDouble value = region[iRow];
164 value *= multiplier[iRow];
165 region[iRow] = value;
166 }
167 }
168 /* Updates one full column (BTRAN) */
updateFullColumnTranspose(CoinIndexedVector & regionSparse) const169 void CoinAbcTypeFactorization::updateFullColumnTranspose(CoinIndexedVector ®ionSparse) const
170 {
171 int numberNonZero = 0;
172 // Should pass in statistics
173 toLongArray(®ionSparse, 0);
174 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
175 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse.getIndices();
176 #ifndef ABC_ORDERED_FACTORIZATION
177 const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
178 for (int iRow = 0; iRow < numberRows_; iRow++) {
179 double value = region[iRow];
180 if (value) {
181 region[iRow] = value * pivotRegion[iRow];
182 regionIndex[numberNonZero++] = iRow;
183 }
184 }
185 regionSparse.setNumElements(numberNonZero);
186 if (!numberNonZero)
187 return;
188 //regionSparse.setNumElements(0);
189 //regionSparse.scan(0,numberRows_);
190 #else
191 permuteInForBtranAndMultiply(regionSparse, true);
192 if (!regionSparse.getNumElements()) {
193 permuteOutForBtran(regionSparse);
194 return;
195 }
196 #endif
197
198 // ******* U
199 // Apply pivot region - could be combined for speed
200 // Can only combine/move inside vector for sparse
201 CoinSimplexInt smallestIndex = pivotLinkedForwardsAddress_[-1];
202 #if ABC_SMALL < 2
203 // copy of code inside transposeU
204 bool goSparse = false;
205 #else
206 #define goSparse false
207 #endif
208 #if ABC_SMALL < 2
209 // Guess at number at end
210 if (gotUCopy()) {
211 assert(btranFullAverageAfterU_);
212 CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * btranFullAverageAfterU_ * twiddleBtranFullFactor1());
213 if (newNumber < sparseThreshold_)
214 goSparse = true;
215 }
216 #endif
217 if (numberNonZero < 40 && (numberNonZero << 4) < numberRows_ && !goSparse) {
218 CoinSimplexInt *COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
219 CoinSimplexInt smallest = numberRowsExtra_;
220 for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
221 CoinSimplexInt iRow = regionIndex[j];
222 if (pivotRowForward[iRow] < smallest) {
223 smallest = pivotRowForward[iRow];
224 smallestIndex = iRow;
225 }
226 }
227 }
228 updateColumnTransposeU(®ionSparse, smallestIndex
229 #if ABC_SMALL < 2
230 ,
231 reinterpret_cast< CoinAbcStatistics & >(btranFullCountInput_)
232 #endif
233 #if ABC_PARALLEL
234 ,
235 0
236 #endif
237 );
238 //row bits here
239 updateColumnTransposeR(®ionSparse
240 #if ABC_SMALL < 2
241 ,
242 reinterpret_cast< CoinAbcStatistics & >(btranFullCountInput_)
243 #endif
244 );
245 // ******* L
246 updateColumnTransposeL(®ionSparse
247 #if ABC_SMALL < 2
248 ,
249 reinterpret_cast< CoinAbcStatistics & >(btranFullCountInput_)
250 #endif
251 #if ABC_PARALLEL
252 ,
253 0
254 #endif
255 );
256 fromLongArray(0);
257 #if ABC_SMALL < 3
258 #ifdef ABC_ORDERED_FACTORIZATION
259 permuteOutForBtran(regionSparse);
260 #endif
261 #if ABC_DEBUG
262 regionSparse.checkClean();
263 #endif
264 #else
265 numberNonZero = 0;
266 for (CoinSimplexInt i = 0; i < numberRows_; i++) {
267 CoinExponent expValue = ABC_EXPONENT(region[i]);
268 if (expValue) {
269 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
270 regionIndex[numberNonZero++] = i;
271 } else {
272 region[i] = 0.0;
273 }
274 }
275 }
276 regionSparse.setNumElements(numberNonZero);
277 #endif
278 }
279 // updateColumnL. Updates part of column (FTRANL)
updateColumnL(CoinIndexedVector * regionSparse,CoinAbcStatistics & statistics,int whichSparse) const280 void CoinAbcTypeFactorization::updateColumnL(CoinIndexedVector *regionSparse
281 #if ABC_SMALL < 2
282 ,
283 CoinAbcStatistics &statistics
284 #endif
285 #if ABC_PARALLEL
286 ,
287 int whichSparse
288 #endif
289 ) const
290 {
291 #if CILK_CONFLICT > 0
292 #if ABC_PARALLEL
293 // for conflicts
294 #if CILK_CONFLICT > 1
295 printf("file %s line %d which %d\n", __FILE__, __LINE__, whichSparse);
296 #endif
297 abc_assert((cilk_conflict & (1 << whichSparse)) == 0);
298 cilk_conflict |= (1 << whichSparse);
299 #else
300 abc_assert((cilk_conflict & (1 << 0)) == 0);
301 cilk_conflict |= (1 << 0);
302 #endif
303 #endif
304 #if ABC_SMALL < 2
305 CoinSimplexInt number = regionSparse->getNumElements();
306 if (factorizationStatistics()) {
307 statistics.numberCounts_++;
308 statistics.countInput_ += number;
309 }
310 #endif
311 if (numberL_) {
312 #if ABC_SMALL < 2
313 int goSparse;
314 // Guess at number at end
315 if (gotSparse()) {
316 double average = statistics.averageAfterL_ * twiddleFactor1S();
317 assert(average);
318 CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(number * average);
319 if (newNumber < sparseThreshold_ && (numberL_ << 2) > newNumber * 1.0 * twiddleFactor2S())
320 goSparse = 1;
321 else if (3 * newNumber < numberRows_)
322 goSparse = 0;
323 else
324 goSparse = -1;
325 } else {
326 goSparse = 0;
327 } //if(goSparse==1) goSparse=0;
328 if (!goSparse) {
329 // densish
330 updateColumnLDensish(regionSparse);
331 } else if (goSparse < 0) {
332 // densish
333 updateColumnLDense(regionSparse);
334 } else {
335 // sparse
336 updateColumnLSparse(regionSparse
337 #if ABC_PARALLEL
338 ,
339 whichSparse
340 #endif
341 );
342 }
343 #else
344 updateColumnLDensish(regionSparse);
345 #endif
346 }
347 #if ABC_SMALL < 4
348 #if ABC_DENSE_CODE > 0
349 if (numberDense_) {
350 instrument_do("CoinAbcFactorizationDense", 0.3 * numberDense_ * numberDense_);
351 //take off list
352 CoinSimplexInt lastSparse = numberRows_ - numberDense_;
353 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
354 CoinFactorizationDouble *COIN_RESTRICT denseArea = denseAreaAddress_;
355 CoinFactorizationDouble *COIN_RESTRICT denseRegion = denseArea + leadingDimension_ * numberDense_;
356 CoinSimplexInt *COIN_RESTRICT densePermute = reinterpret_cast< CoinSimplexInt * >(denseRegion + FACTOR_CPU * numberDense_);
357 //for (int i=0;i<numberDense_;i++) {
358 //printf("perm %d %d\n",i,densePermute[i]);
359 //assert (densePermute[i]>=0&&densePermute[i]<numberRows_);
360 //}
361 #if ABC_PARALLEL
362 if (whichSparse)
363 denseRegion += whichSparse * numberDense_;
364 //printf("PP %d %d %s region %x\n",whichSparse,__LINE__,__FILE__,denseRegion);
365 #endif
366 CoinFactorizationDouble *COIN_RESTRICT densePut = denseRegion - lastSparse;
367 CoinZeroN(denseRegion, numberDense_);
368 bool doDense = false;
369 #if ABC_SMALL < 3
370 #if ABC_DENSE_CODE == 2
371 //short * COIN_RESTRICT forFtran = reinterpret_cast<short *>(densePermute+numberDense_)-lastSparse;
372 short *COIN_RESTRICT forFtran2 = reinterpret_cast< short * >(densePermute + numberDense_) + 2 * numberDense_;
373 #else
374 const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
375 #endif
376 CoinSimplexInt number = regionSparse->getNumElements();
377 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
378 CoinSimplexInt i = 0;
379 while (i < number) {
380 CoinSimplexInt iRow = regionIndex[i];
381 #if ABC_DENSE_CODE == 2
382 int kRow = forFtran2[iRow];
383 if (kRow != -1) {
384 doDense = true;
385 regionIndex[i] = regionIndex[--number];
386 denseRegion[kRow] = region[iRow];
387 #else
388 CoinSimplexInt jRow = pivotLBackwardOrder[iRow];
389 if (jRow >= lastSparse) {
390 doDense = true;
391 regionIndex[i] = regionIndex[--number];
392 densePut[jRow] = region[iRow];
393 #endif
394 region[iRow] = 0.0;
395 } else {
396 i++;
397 }
398 }
399 #else
400 CoinSimplexInt *COIN_RESTRICT forFtran = densePermute + numberDense_ - lastSparse;
401 const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
402 for (CoinSimplexInt jRow = lastSparse; jRow < numberRows_; jRow++) {
403 CoinSimplexInt iRow = pivotLOrder[jRow];
404 if (region[iRow]) {
405 doDense = true;
406 //densePut[jRow]=region[iRow];
407 #if ABC_DENSE_CODE == 2
408 int kRow = forFtran[jRow];
409 denseRegion[kRow] = region[iRow];
410 #endif
411 region[iRow] = 0.0;
412 }
413 }
414 #endif
415 if (doDense) {
416 #if ABC_DENSE_CODE == 1
417 #ifndef CLAPACK
418 char trans = 'N';
419 CoinSimplexInt ione = 1;
420 CoinSimplexInt info;
421 F77_FUNC(dgetrs, DGETRS)
422 (&trans, &numberDense_, &ione, denseArea, &leadingDimension_,
423 densePermute, denseRegion, &numberDense_, &info, 1);
424 #else
425 clapack_dgetrs(CblasColMajor, CblasNoTrans, numberDense_, 1,
426 denseArea, leadingDimension_, densePermute, denseRegion, numberDense_);
427 #endif
428 #elif ABC_DENSE_CODE == 2
429 CoinAbcDgetrs('N', numberDense_, denseArea, denseRegion);
430 #endif
431 const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
432 for (CoinSimplexInt i = lastSparse; i < numberRows_; i++) {
433 if (!TEST_LESS_THAN_TOLERANCE(densePut[i])) {
434 #ifndef ABC_ORDERED_FACTORIZATION
435 CoinSimplexInt iRow = pivotLOrder[i];
436 #else
437 CoinSimplexInt iRow = i;
438 #endif
439 region[iRow] = densePut[i];
440 #if ABC_SMALL < 3
441 regionIndex[number++] = iRow;
442 #endif
443 }
444 }
445 #if ABC_SMALL < 3
446 regionSparse->setNumElements(number);
447 #endif
448 }
449 }
450 //printRegion(*regionSparse,"After FtranL");
451 #endif
452 #endif
453 #if ABC_SMALL < 2
454 if (factorizationStatistics())
455 statistics.countAfterL_ += regionSparse->getNumElements();
456 #endif
457 #if CILK_CONFLICT > 0
458 #if ABC_PARALLEL
459 // for conflicts
460 abc_assert((cilk_conflict & (1 << whichSparse)) != 0);
461 cilk_conflict &= ~(1 << whichSparse);
462 #else
463 abc_assert((cilk_conflict & (1 << 0)) != 0);
464 cilk_conflict &= ~(1 << 0);
465 #endif
466 #endif
467 }
468 #define UNROLL 2
469 #define INLINE_IT
470 //#define INLINE_IT2
471 inline void scatterUpdateInline(CoinSimplexInt number,
472 CoinFactorizationDouble pivotValue,
473 const CoinFactorizationDouble *COIN_RESTRICT thisElement,
474 const CoinSimplexInt *COIN_RESTRICT thisIndex,
475 CoinFactorizationDouble *COIN_RESTRICT region)
476 {
477 #if UNROLL == 0
478 for (CoinBigIndex j = number - 1; j >= 0; j--) {
479 CoinSimplexInt iRow = thisIndex[j];
480 CoinFactorizationDouble regionValue = region[iRow];
481 CoinFactorizationDouble value = thisElement[j];
482 assert(value);
483 region[iRow] = regionValue - value * pivotValue;
484 }
485 #elif UNROLL == 1
486 if ((number & 1) != 0) {
487 number--;
488 CoinSimplexInt iRow = thisIndex[number];
489 CoinFactorizationDouble regionValue = region[iRow];
490 CoinFactorizationDouble value = thisElement[number];
491 region[iRow] = regionValue - value * pivotValue;
492 }
493 for (CoinBigIndex j = number - 1; j >= 0; j -= 2) {
494 CoinSimplexInt iRow0 = thisIndex[j];
495 CoinSimplexInt iRow1 = thisIndex[j - 1];
496 CoinFactorizationDouble regionValue0 = region[iRow0];
497 CoinFactorizationDouble regionValue1 = region[iRow1];
498 region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
499 region[iRow1] = regionValue1 - thisElement[j - 1] * pivotValue;
500 }
501 #elif UNROLL == 2
502 CoinSimplexInt iRow0;
503 CoinSimplexInt iRow1;
504 CoinFactorizationDouble regionValue0;
505 CoinFactorizationDouble regionValue1;
506 switch (static_cast< unsigned int >(number)) {
507 case 0:
508 break;
509 case 1:
510 iRow0 = thisIndex[0];
511 regionValue0 = region[iRow0];
512 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
513 break;
514 case 2:
515 iRow0 = thisIndex[0];
516 iRow1 = thisIndex[1];
517 regionValue0 = region[iRow0];
518 regionValue1 = region[iRow1];
519 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
520 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
521 break;
522 case 3:
523 iRow0 = thisIndex[0];
524 iRow1 = thisIndex[1];
525 regionValue0 = region[iRow0];
526 regionValue1 = region[iRow1];
527 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
528 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
529 iRow0 = thisIndex[2];
530 regionValue0 = region[iRow0];
531 region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
532 break;
533 case 4:
534 iRow0 = thisIndex[0];
535 iRow1 = thisIndex[1];
536 regionValue0 = region[iRow0];
537 regionValue1 = region[iRow1];
538 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
539 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
540 iRow0 = thisIndex[2];
541 iRow1 = thisIndex[3];
542 regionValue0 = region[iRow0];
543 regionValue1 = region[iRow1];
544 region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
545 region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
546 break;
547 case 5:
548 iRow0 = thisIndex[0];
549 iRow1 = thisIndex[1];
550 regionValue0 = region[iRow0];
551 regionValue1 = region[iRow1];
552 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
553 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
554 iRow0 = thisIndex[2];
555 iRow1 = thisIndex[3];
556 regionValue0 = region[iRow0];
557 regionValue1 = region[iRow1];
558 region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
559 region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
560 iRow0 = thisIndex[4];
561 regionValue0 = region[iRow0];
562 region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
563 break;
564 case 6:
565 iRow0 = thisIndex[0];
566 iRow1 = thisIndex[1];
567 regionValue0 = region[iRow0];
568 regionValue1 = region[iRow1];
569 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
570 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
571 iRow0 = thisIndex[2];
572 iRow1 = thisIndex[3];
573 regionValue0 = region[iRow0];
574 regionValue1 = region[iRow1];
575 region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
576 region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
577 iRow0 = thisIndex[4];
578 iRow1 = thisIndex[5];
579 regionValue0 = region[iRow0];
580 regionValue1 = region[iRow1];
581 region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
582 region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
583 break;
584 case 7:
585 iRow0 = thisIndex[0];
586 iRow1 = thisIndex[1];
587 regionValue0 = region[iRow0];
588 regionValue1 = region[iRow1];
589 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
590 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
591 iRow0 = thisIndex[2];
592 iRow1 = thisIndex[3];
593 regionValue0 = region[iRow0];
594 regionValue1 = region[iRow1];
595 region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
596 region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
597 iRow0 = thisIndex[4];
598 iRow1 = thisIndex[5];
599 regionValue0 = region[iRow0];
600 regionValue1 = region[iRow1];
601 region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
602 region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
603 iRow0 = thisIndex[6];
604 regionValue0 = region[iRow0];
605 region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
606 break;
607 case 8:
608 iRow0 = thisIndex[0];
609 iRow1 = thisIndex[1];
610 regionValue0 = region[iRow0];
611 regionValue1 = region[iRow1];
612 region[iRow0] = regionValue0 - thisElement[0] * pivotValue;
613 region[iRow1] = regionValue1 - thisElement[1] * pivotValue;
614 iRow0 = thisIndex[2];
615 iRow1 = thisIndex[3];
616 regionValue0 = region[iRow0];
617 regionValue1 = region[iRow1];
618 region[iRow0] = regionValue0 - thisElement[2] * pivotValue;
619 region[iRow1] = regionValue1 - thisElement[3] * pivotValue;
620 iRow0 = thisIndex[4];
621 iRow1 = thisIndex[5];
622 regionValue0 = region[iRow0];
623 regionValue1 = region[iRow1];
624 region[iRow0] = regionValue0 - thisElement[4] * pivotValue;
625 region[iRow1] = regionValue1 - thisElement[5] * pivotValue;
626 iRow0 = thisIndex[6];
627 iRow1 = thisIndex[7];
628 regionValue0 = region[iRow0];
629 regionValue1 = region[iRow1];
630 region[iRow0] = regionValue0 - thisElement[6] * pivotValue;
631 region[iRow1] = regionValue1 - thisElement[7] * pivotValue;
632 break;
633 default:
634 if ((number & 1) != 0) {
635 number--;
636 CoinSimplexInt iRow = thisIndex[number];
637 CoinFactorizationDouble regionValue = region[iRow];
638 CoinFactorizationDouble value = thisElement[number];
639 region[iRow] = regionValue - value * pivotValue;
640 }
641 for (CoinBigIndex j = number - 1; j >= 0; j -= 2) {
642 CoinSimplexInt iRow0 = thisIndex[j];
643 CoinSimplexInt iRow1 = thisIndex[j - 1];
644 CoinFactorizationDouble regionValue0 = region[iRow0];
645 CoinFactorizationDouble regionValue1 = region[iRow1];
646 region[iRow0] = regionValue0 - thisElement[j] * pivotValue;
647 region[iRow1] = regionValue1 - thisElement[j - 1] * pivotValue;
648 }
649 break;
650 }
651 #endif
652 }
653 inline CoinFactorizationDouble gatherUpdate(CoinSimplexInt number,
654 const CoinFactorizationDouble *COIN_RESTRICT thisElement,
655 const CoinSimplexInt *COIN_RESTRICT thisIndex,
656 CoinFactorizationDouble *COIN_RESTRICT region)
657 {
658 CoinFactorizationDouble pivotValue = 0.0;
659 for (CoinBigIndex j = 0; j < number; j++) {
660 CoinFactorizationDouble value = thisElement[j];
661 CoinSimplexInt jRow = thisIndex[j];
662 value *= region[jRow];
663 pivotValue -= value;
664 }
665 return pivotValue;
666 }
667 // Updates part of column (FTRANL) when densish
668 void CoinAbcTypeFactorization::updateColumnLDensish(CoinIndexedVector *regionSparse) const
669 {
670 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
671 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
672 CoinSimplexInt number = regionSparse->getNumElements();
673 #if ABC_SMALL < 3
674 CoinSimplexInt numberNonZero = 0;
675 #endif
676
677 const CoinBigIndex *COIN_RESTRICT startColumn = startColumnLAddress_;
678 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowLAddress_;
679 const CoinFactorizationDouble *COIN_RESTRICT element = elementLAddress_;
680 const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
681 const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
682 CoinSimplexInt last = numberRows_;
683 assert(last == baseL_ + numberL_);
684 #if ABC_DENSE_CODE > 0
685 //can take out last bit of sparse L as empty
686 last -= numberDense_;
687 #endif
688 CoinSimplexInt smallestIndex = numberRowsExtra_;
689 // do easy ones
690 for (CoinSimplexInt k = 0; k < number; k++) {
691 CoinSimplexInt iPivot = regionIndex[k];
692 #ifndef ABC_ORDERED_FACTORIZATION
693 CoinSimplexInt jPivot = pivotLBackwardOrder[iPivot];
694 #else
695 CoinSimplexInt jPivot = iPivot;
696 #endif
697 #if ABC_SMALL < 3
698 if (jPivot >= baseL_)
699 smallestIndex = CoinMin(jPivot, smallestIndex);
700 else
701 regionIndex[numberNonZero++] = iPivot;
702 #else
703 if (jPivot >= baseL_)
704 smallestIndex = CoinMin(jPivot, smallestIndex);
705 #endif
706 }
707 instrument_start("CoinAbcFactorizationUpdateLDensish", number + (last - smallestIndex));
708 // now others
709 for (CoinSimplexInt k = smallestIndex; k < last; k++) {
710 #if 0
711 for (int j=0;j<numberRows_;j+=10) {
712 for (int jj=j;jj<CoinMin(j+10,numberRows_);jj++)
713 printf("%g ",region[jj]);
714 printf("\n");
715 }
716 #endif
717 #ifndef ABC_ORDERED_FACTORIZATION
718 CoinSimplexInt i = pivotLOrder[k];
719 #else
720 CoinSimplexInt i = k;
721 #endif
722 CoinExponent expValue = ABC_EXPONENT(region[i]);
723 if (expValue) {
724 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
725 CoinBigIndex start = startColumn[k];
726 CoinBigIndex end = startColumn[k + 1];
727 instrument_add(end - start);
728 if (TEST_INT_NONZERO(end - start)) {
729 CoinFactorizationDouble pivotValue = region[i];
730 #ifndef INLINE_IT
731 for (CoinBigIndex j = start; j < end; j++) {
732 CoinSimplexInt iRow = indexRow[j];
733 CoinFactorizationDouble result = region[iRow];
734 CoinFactorizationDouble value = element[j];
735 region[iRow] = result - value * pivotValue;
736 }
737 #else
738 CoinAbcScatterUpdate(end - start, pivotValue, element + start, indexRow + start, region);
739 #endif
740 }
741 #if ABC_SMALL < 3
742 regionIndex[numberNonZero++] = i;
743 } else {
744 region[i] = 0.0;
745 #endif
746 }
747 }
748 }
749 // and dense
750 #if ABC_SMALL < 3
751 for (CoinSimplexInt k = last; k < numberRows_; k++) {
752 #ifndef ABC_ORDERED_FACTORIZATION
753 CoinSimplexInt i = pivotLOrder[k];
754 #else
755 CoinSimplexInt i = k;
756 #endif
757 CoinExponent expValue = ABC_EXPONENT(region[i]);
758 if (expValue) {
759 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
760 regionIndex[numberNonZero++] = i;
761 } else {
762 region[i] = 0.0;
763 }
764 }
765 }
766 regionSparse->setNumElements(numberNonZero);
767 #endif
768 instrument_end();
769 }
770 // Updates part of column (FTRANL) when dense
771 void CoinAbcTypeFactorization::updateColumnLDense(CoinIndexedVector *regionSparse) const
772 {
773 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
774 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
775 CoinSimplexInt number = regionSparse->getNumElements();
776 #if ABC_SMALL < 3
777 CoinSimplexInt numberNonZero = 0;
778 #endif
779
780 const CoinBigIndex *COIN_RESTRICT startColumn = startColumnLAddress_;
781 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowLAddress_;
782 const CoinFactorizationDouble *COIN_RESTRICT element = elementLAddress_;
783 const CoinSimplexInt *COIN_RESTRICT pivotLOrder = pivotLOrderAddress_;
784 const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
785 CoinSimplexInt last = numberRows_;
786 assert(last == baseL_ + numberL_);
787 #if ABC_DENSE_CODE > 0
788 //can take out last bit of sparse L as empty
789 last -= numberDense_;
790 #endif
791 CoinSimplexInt smallestIndex = numberRowsExtra_;
792 // do easy ones
793 for (CoinSimplexInt k = 0; k < number; k++) {
794 CoinSimplexInt iPivot = regionIndex[k];
795 #ifndef ABC_ORDERED_FACTORIZATION
796 CoinSimplexInt jPivot = pivotLBackwardOrder[iPivot];
797 #else
798 CoinSimplexInt jPivot = iPivot;
799 #endif
800 #if ABC_SMALL < 3
801 if (jPivot >= baseL_)
802 smallestIndex = CoinMin(jPivot, smallestIndex);
803 else
804 regionIndex[numberNonZero++] = iPivot;
805 #else
806 if (jPivot >= baseL_)
807 smallestIndex = CoinMin(jPivot, smallestIndex);
808 #endif
809 }
810 instrument_start("CoinAbcFactorizationUpdateLDensish", number + (last - smallestIndex));
811 // now others
812 for (CoinSimplexInt k = smallestIndex; k < last; k++) {
813 #ifndef ABC_ORDERED_FACTORIZATION
814 CoinSimplexInt i = pivotLOrder[k];
815 #else
816 CoinSimplexInt i = k;
817 #endif
818 CoinExponent expValue = ABC_EXPONENT(region[i]);
819 if (expValue) {
820 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
821 CoinBigIndex start = startColumn[k];
822 CoinBigIndex end = startColumn[k + 1];
823 instrument_add(end - start);
824 if (TEST_INT_NONZERO(end - start)) {
825 CoinFactorizationDouble pivotValue = region[i];
826 #ifndef INLINE_IT
827 for (CoinBigIndex j = start; j < end; j++) {
828 CoinSimplexInt iRow = indexRow[j];
829 CoinFactorizationDouble result = region[iRow];
830 CoinFactorizationDouble value = element[j];
831 region[iRow] = result - value * pivotValue;
832 }
833 #else
834 CoinAbcScatterUpdate(end - start, pivotValue, element + start, indexRow + start, region);
835 #endif
836 }
837 #if ABC_SMALL < 3
838 regionIndex[numberNonZero++] = i;
839 } else {
840 region[i] = 0.0;
841 #endif
842 }
843 }
844 }
845 // and dense
846 #if ABC_SMALL < 3
847 for (CoinSimplexInt k = last; k < numberRows_; k++) {
848 #ifndef ABC_ORDERED_FACTORIZATION
849 CoinSimplexInt i = pivotLOrder[k];
850 #else
851 CoinSimplexInt i = k;
852 #endif
853 CoinExponent expValue = ABC_EXPONENT(region[i]);
854 if (expValue) {
855 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
856 regionIndex[numberNonZero++] = i;
857 } else {
858 region[i] = 0.0;
859 }
860 }
861 }
862 regionSparse->setNumElements(numberNonZero);
863 #endif
864 instrument_end();
865 }
866 #if ABC_SMALL < 2
867 // Updates part of column (FTRANL) when sparse
868 void CoinAbcTypeFactorization::updateColumnLSparse(CoinIndexedVector *regionSparse
869 #if ABC_PARALLEL
870 ,
871 int whichSparse
872 #endif
873 ) const
874 {
875 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
876 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
877 CoinSimplexInt number = regionSparse->getNumElements();
878 CoinSimplexInt numberNonZero;
879
880 numberNonZero = 0;
881
882 const CoinBigIndex *COIN_RESTRICT startColumn = startColumnLAddress_;
883 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowLAddress_;
884 const CoinFactorizationDouble *COIN_RESTRICT element = elementLAddress_;
885 const CoinSimplexInt *COIN_RESTRICT pivotLBackwardOrder = permuteAddress_;
886 CoinSimplexInt nList;
887 // use sparse_ as temporary area
888 // need to redo if this fails (just means using CoinAbcStack to compute sizes)
889 assert(sizeof(CoinSimplexInt) == sizeof(CoinBigIndex));
890 CoinAbcStack *COIN_RESTRICT stackList = reinterpret_cast< CoinAbcStack * >(sparseAddress_);
891 CoinSimplexInt *COIN_RESTRICT list = listAddress_;
892 #define DO_CHAR1 1
893 #if DO_CHAR1 // CHAR
894 CoinCheckZero *COIN_RESTRICT mark = markRowAddress_;
895 #else
896 //BIT
897 CoinSimplexUnsignedInt *COIN_RESTRICT mark = reinterpret_cast< CoinSimplexUnsignedInt * >(markRowAddress_);
898 #endif
899 #if ABC_PARALLEL
900 //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
901 if (whichSparse) {
902 //printf("alternative sparse\n");
903 int addAmount = whichSparse * sizeSparseArray_;
904 stackList = reinterpret_cast< CoinAbcStack * >(reinterpret_cast< char * >(stackList) + addAmount);
905 list = reinterpret_cast< CoinSimplexInt * >(reinterpret_cast< char * >(list) + addAmount);
906 #if DO_CHAR1 // CHAR
907 mark = reinterpret_cast< CoinCheckZero * >(reinterpret_cast< char * >(mark) + addAmount);
908 #else
909 mark = reinterpret_cast< CoinSimplexUnsignedInt * >(reinterpret_cast< char * >(mark) + addAmount);
910 #endif
911 }
912 #endif
913 // mark known to be zero
914 #ifdef COIN_DEBUG
915 #if DO_CHAR1 // CHAR
916 for (CoinSimplexInt i = 0; i < maximumRowsExtra_; i++) {
917 assert(!mark[i]);
918 }
919 #else
920 //BIT
921 for (CoinSimplexInt i = 0; i< numberRows_ >> COINFACTORIZATION_SHIFT_PER_INT; i++) {
922 assert(!mark[i]);
923 }
924 #endif
925 #endif
926 nList = 0;
927 for (CoinSimplexInt k = 0; k < number; k++) {
928 CoinSimplexInt kPivot = regionIndex[k];
929 #ifndef ABC_ORDERED_FACTORIZATION
930 CoinSimplexInt iPivot = pivotLBackwardOrder[kPivot];
931 #else
932 CoinSimplexInt iPivot = kPivot;
933 #endif
934 if (iPivot >= baseL_) {
935 #if DO_CHAR1 // CHAR
936 CoinCheckZero mark_k = mark[kPivot];
937 #else
938 CoinSimplexUnsignedInt wordK = kPivot >> COINFACTORIZATION_SHIFT_PER_INT;
939 CoinSimplexUnsignedInt bitK = (kPivot & COINFACTORIZATION_MASK_PER_INT);
940 CoinSimplexUnsignedInt mark_k = (mark[wordK] >> bitK) & 1;
941 #endif
942 if (!mark_k) {
943 CoinBigIndex j = startColumn[iPivot + 1] - 1;
944 stackList[0].stack = kPivot;
945 CoinBigIndex start = startColumn[iPivot];
946 stackList[0].next = j;
947 stackList[0].start = start;
948 CoinSimplexInt nStack = 0;
949 while (nStack >= 0) {
950 /* take off stack */
951 #ifndef ABC_ORDERED_FACTORIZATION
952 iPivot = pivotLBackwardOrder[kPivot]; // put startColumn on stackstuff
953 #else
954 iPivot = kPivot;
955 #endif
956 #if 1 //0 // what went wrong??
957 CoinBigIndex startThis = startColumn[iPivot];
958 for (; j >= startThis; j--) {
959 CoinSimplexInt jPivot = indexRow[j];
960 #if DO_CHAR1 // CHAR
961 CoinCheckZero mark_j = mark[jPivot];
962 #else
963 CoinSimplexUnsignedInt word0 = jPivot >> COINFACTORIZATION_SHIFT_PER_INT;
964 CoinSimplexUnsignedInt bit0 = (jPivot & COINFACTORIZATION_MASK_PER_INT);
965 CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 1;
966 #endif
967 if (!mark_j) {
968 #if DO_CHAR1 // CHAR
969 mark[jPivot] = 1;
970 #else
971 mark[word0] |= (1 << bit0);
972 #endif
973 /* and new one */
974 kPivot = jPivot;
975 break;
976 }
977 }
978 if (j >= startThis) {
979 /* put back on stack */
980 stackList[nStack].next = j - 1;
981 #ifndef ABC_ORDERED_FACTORIZATION
982 iPivot = pivotLBackwardOrder[kPivot];
983 #else
984 iPivot = kPivot;
985 #endif
986 j = startColumn[iPivot + 1] - 1;
987 nStack++;
988 stackList[nStack].stack = kPivot;
989 assert(kPivot < numberRowsExtra_);
990 CoinBigIndex start = startColumn[iPivot];
991 stackList[nStack].next = j;
992 stackList[nStack].start = start;
993 #else
994 if (j >= startColumn[iPivot] /*stackList[nStack].start*/) {
995 CoinSimplexInt jPivot = indexRow[j--];
996 /* put back on stack */
997 stackList[nStack].next = j;
998 #if DO_CHAR1 // CHAR
999 CoinCheckZero mark_j = mark[jPivot];
1000 #else
1001 CoinSimplexUnsignedInt word0 = jPivot >> COINFACTORIZATION_SHIFT_PER_INT;
1002 CoinSimplexUnsignedInt bit0 = (jPivot & COINFACTORIZATION_MASK_PER_INT);
1003 CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 1;
1004 #endif
1005 if (!mark_j) {
1006 /* and new one */
1007 kPivot = jPivot;
1008 #ifndef ABC_ORDERED_FACTORIZATION
1009 iPivot = pivotLBackwardOrder[kPivot];
1010 #else
1011 iPivot = kPivot;
1012 #endif
1013 j = startColumn[iPivot + 1] - 1;
1014 nStack++;
1015 stackList[nStack].stack = kPivot;
1016 assert(kPivot < numberRowsExtra_);
1017 CoinBigIndex start = startColumn[iPivot];
1018 stackList[nStack].next = j;
1019 stackList[nStack].start = start;
1020 #if DO_CHAR1 // CHAR
1021 mark[jPivot] = 1;
1022 #else
1023 mark[word0] |= (1 << bit0);
1024 #endif
1025 }
1026 #endif
1027 } else {
1028 /* finished so mark */
1029 list[nList++] = kPivot;
1030 #if DO_CHAR1 // CHAR
1031 mark[kPivot] = 1;
1032 #else
1033 CoinSimplexUnsignedInt wordB = kPivot >> COINFACTORIZATION_SHIFT_PER_INT;
1034 CoinSimplexUnsignedInt bitB = (kPivot & COINFACTORIZATION_MASK_PER_INT);
1035 mark[wordB] |= (1 << bitB);
1036 #endif
1037 --nStack;
1038 if (nStack >= 0) {
1039 kPivot = stackList[nStack].stack;
1040 #ifndef ABC_ORDERED_FACTORIZATION
1041 iPivot = pivotLBackwardOrder[kPivot];
1042 #else
1043 iPivot = kPivot;
1044 #endif
1045 assert(kPivot < numberRowsExtra_);
1046 j = stackList[nStack].next;
1047 }
1048 }
1049 }
1050 }
1051 } else {
1052 if (!TEST_LESS_THAN_TOLERANCE(region[kPivot])) {
1053 // just put on list
1054 regionIndex[numberNonZero++] = kPivot;
1055 } else {
1056 region[kPivot] = 0.0;
1057 }
1058 }
1059 }
1060 #if 0
1061 CoinSimplexDouble temp[20000];
1062 {
1063 memcpy(temp,region,numberRows_*sizeof(CoinSimplexDouble));
1064 for (CoinSimplexInt k = 0; k < numberRows_; k++ ) {
1065 CoinSimplexInt i=pivotLOrder[k];
1066 CoinFactorizationDouble pivotValue = temp[i];
1067
1068 if ( !TEST_LESS_THAN_TOLERANCE(pivotValue) ) {
1069 CoinBigIndex start = startColumn[k];
1070 CoinBigIndex end = startColumn[k + 1];
1071 for (CoinBigIndex j = start; j < end; j ++ ) {
1072 CoinSimplexInt iRow = indexRow[j];
1073 CoinFactorizationDouble result = temp[iRow];
1074 CoinFactorizationDouble value = element[j];
1075
1076 temp[iRow] = result - value * pivotValue;
1077 }
1078 } else {
1079 temp[i] = 0.0;
1080 }
1081 }
1082 }
1083 #endif
1084 instrument_start("CoinAbcFactorizationUpdateLSparse", number + nList);
1085 #if 0 //ndef NDEBUG
1086 char * test = new char[numberRows_];
1087 memset(test,0,numberRows_);
1088 for (int i=0;i<numberRows_;i++) {
1089 if (region[i])
1090 test[i]=-1;
1091 }
1092 for (CoinSimplexInt i=nList-1;i>=0;i--) {
1093 CoinSimplexInt iPivot = list[i];
1094 test[iPivot]=1;
1095 CoinSimplexInt kPivot=pivotLBackwardOrder[iPivot];
1096 CoinBigIndex start=startColumn[kPivot];
1097 CoinBigIndex end=startColumn[kPivot+1];
1098 for (CoinBigIndex j = start;j < end; j ++ ) {
1099 CoinSimplexInt iRow = indexRow[j];
1100 assert (test[iRow]<=0);
1101 }
1102 }
1103 delete [] test;
1104 #endif
1105 for (CoinSimplexInt i = nList - 1; i >= 0; i--) {
1106 CoinSimplexInt iPivot = list[i];
1107 #ifndef ABC_ORDERED_FACTORIZATION
1108 CoinSimplexInt kPivot = pivotLBackwardOrder[iPivot];
1109 #else
1110 CoinSimplexInt kPivot = iPivot;
1111 #endif
1112 #if DO_CHAR1 // CHAR
1113 mark[iPivot] = 0;
1114 #else
1115 CoinSimplexUnsignedInt word0 = iPivot >> COINFACTORIZATION_SHIFT_PER_INT;
1116 //CoinSimplexUnsignedInt bit0 = (iPivot & COINFACTORIZATION_MASK_PER_INT);
1117 //mark[word0]&=(~(1<<bit0));
1118 mark[word0] = 0;
1119 #endif
1120 if (!TEST_LESS_THAN_TOLERANCE(region[iPivot])) {
1121 CoinFactorizationDouble pivotValue = region[iPivot];
1122 regionIndex[numberNonZero++] = iPivot;
1123 CoinBigIndex start = startColumn[kPivot];
1124 CoinBigIndex end = startColumn[kPivot + 1];
1125 instrument_add(end - start);
1126 if (TEST_INT_NONZERO(end - start)) {
1127 #ifndef INLINE_IT
1128 for (CoinBigIndex j = start; j < end; j++) {
1129 CoinSimplexInt iRow = indexRow[j];
1130 CoinFactorizationDouble value = element[j];
1131 region[iRow] -= value * pivotValue;
1132 }
1133 #else
1134 CoinAbcScatterUpdate(end - start, pivotValue, element + start, indexRow + start, region);
1135 #endif
1136 }
1137 } else {
1138 region[iPivot] = 0.0;
1139 }
1140 }
1141 #if 0
1142 {
1143 for (CoinSimplexInt k = 0; k < numberRows_; k++ ) {
1144 assert (fabs(region[k]-temp[k])<1.0e-1);
1145 }
1146 }
1147 #endif
1148 regionSparse->setNumElements(numberNonZero);
1149 instrument_end_and_adjust(1.3);
1150 }
1151 #endif
1152 #if FACTORIZATION_STATISTICS
1153 extern double twoThresholdX;
1154 #endif
1155 CoinSimplexInt
1156 CoinAbcTypeFactorization::updateTwoColumnsFT(CoinIndexedVector ®ionFTX,
1157 CoinIndexedVector ®ionOtherX)
1158 {
1159 CoinIndexedVector *regionFT = ®ionFTX;
1160 CoinIndexedVector *regionOther = ®ionOtherX;
1161 #if ABC_DEBUG
1162 regionFT->checkClean();
1163 regionOther->checkClean();
1164 #endif
1165 toLongArray(regionFT, 0);
1166 toLongArray(regionOther, 1);
1167 #ifdef ABC_ORDERED_FACTORIZATION
1168 // Permute in for Ftran
1169 permuteInForFtran(regionFTX);
1170 permuteInForFtran(regionOtherX);
1171 #endif
1172 CoinSimplexInt numberNonZero = regionOther->getNumElements();
1173 CoinSimplexInt numberNonZeroFT = regionFT->getNumElements();
1174 // ******* L
1175 updateColumnL(regionFT
1176 #if ABC_SMALL < 2
1177 ,
1178 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
1179 #endif
1180 );
1181 updateColumnL(regionOther
1182 #if ABC_SMALL < 2
1183 ,
1184 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
1185 #endif
1186 );
1187 //row bits here
1188 updateColumnR(regionFT
1189 #if ABC_SMALL < 2
1190 ,
1191 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
1192 #endif
1193 );
1194 updateColumnR(regionOther
1195 #if ABC_SMALL < 2
1196 ,
1197 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
1198 #endif
1199 );
1200 bool doFT = storeFT(regionFT);
1201 //update counts
1202 // ******* U - see if densish
1203 #if ABC_SMALL < 2
1204 // Guess at number at end
1205 CoinSimplexInt goSparse = 0;
1206 if (gotSparse()) {
1207 CoinSimplexInt numberNonZero = (regionOther->getNumElements() + regionFT->getNumElements()) >> 1;
1208 double average = 0.25 * (ftranAverageAfterL_ * twiddleFtranFactor1() + ftranFTAverageAfterL_ * twiddleFtranFTFactor1());
1209 assert(average);
1210 CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * average);
1211 if (newNumber < sparseThreshold_)
1212 goSparse = 2;
1213 }
1214 #if FACTORIZATION_STATISTICS
1215 int twoThreshold = twoThresholdX;
1216 #else
1217 #define twoThreshold 1000
1218 #endif
1219 #else
1220 #define goSparse false
1221 #define twoThreshold 1000
1222 #endif
1223 if (!goSparse && numberRows_ < twoThreshold) {
1224 CoinFactorizationDouble *COIN_RESTRICT arrayFT = denseVector(regionFT);
1225 CoinSimplexInt *COIN_RESTRICT indexFT = regionFT->getIndices();
1226 CoinFactorizationDouble *COIN_RESTRICT arrayUpdate = denseVector(regionOther);
1227 CoinSimplexInt *COIN_RESTRICT indexUpdate = regionOther->getIndices();
1228 updateTwoColumnsUDensish(numberNonZeroFT, arrayFT, indexFT,
1229 numberNonZero, arrayUpdate, indexUpdate);
1230 regionFT->setNumElements(numberNonZeroFT);
1231 regionOther->setNumElements(numberNonZero);
1232 } else {
1233 // sparse
1234 updateColumnU(regionFT
1235 #if ABC_SMALL < 2
1236 ,
1237 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
1238 #endif
1239 );
1240 numberNonZeroFT = regionFT->getNumElements();
1241 updateColumnU(regionOther
1242 #if ABC_SMALL < 2
1243 ,
1244 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
1245 #endif
1246 );
1247 numberNonZero = regionOther->getNumElements();
1248 }
1249 fromLongArray(0);
1250 fromLongArray(1);
1251 #if ABC_DEBUG
1252 regionFT->checkClean();
1253 regionOther->checkClean();
1254 #endif
1255 if (doFT)
1256 return numberNonZeroFT;
1257 else
1258 return -numberNonZeroFT;
1259 }
1260 // Updates part of 2 columns (FTRANU) real work
1261 void CoinAbcTypeFactorization::updateTwoColumnsUDensish(
1262 CoinSimplexInt &numberNonZero1,
1263 CoinFactorizationDouble *COIN_RESTRICT region1,
1264 CoinSimplexInt *COIN_RESTRICT index1,
1265 CoinSimplexInt &numberNonZero2,
1266 CoinFactorizationDouble *COIN_RESTRICT region2,
1267 CoinSimplexInt *COIN_RESTRICT index2) const
1268 {
1269 #ifdef ABC_USE_FUNCTION_POINTERS
1270 scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1271 CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1272 #else
1273 const CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1274 const CoinBigIndex *COIN_RESTRICT startColumn = startColumnUAddress_;
1275 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1276 const CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1277 #endif
1278 CoinSimplexInt numberNonZeroA = 0;
1279 CoinSimplexInt numberNonZeroB = 0;
1280 const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1281 const CoinSimplexInt *COIN_RESTRICT pivotLinked = pivotLinkedBackwardsAddress_;
1282 CoinSimplexInt jRow = pivotLinked[numberRows_];
1283 instrument_start("CoinAbcFactorizationUpdateTwoUDense", 2 * numberRows_);
1284 #if 1 //def DONT_USE_SLACKS
1285 while (jRow != -1 /*lastSlack_*/) {
1286 #else
1287 // would need extra code
1288 while (jRow != lastSlack_) {
1289 #endif
1290 bool nonZero2 = (TEST_DOUBLE_NONZERO(region2[jRow]));
1291 bool nonZero1 = (TEST_DOUBLE_NONZERO(region1[jRow]));
1292 #ifndef ABC_USE_FUNCTION_POINTERS
1293 int numberIn = numberInColumn[jRow];
1294 #else
1295 scatterStruct &COIN_RESTRICT scatter = scatterPointer[jRow];
1296 CoinSimplexInt numberIn = scatter.number;
1297 #endif
1298 if (nonZero1 || nonZero2) {
1299 #ifndef ABC_USE_FUNCTION_POINTERS
1300 CoinBigIndex start = startColumn[jRow];
1301 const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
1302 const CoinSimplexInt *COIN_RESTRICT thisIndex = indexRow + start;
1303 #else
1304 const CoinFactorizationDouble *COIN_RESTRICT thisElement = elementUColumnPlus + scatter.offset;
1305 const CoinSimplexInt *COIN_RESTRICT thisIndex = reinterpret_cast< const CoinSimplexInt * >(thisElement + numberIn);
1306 #endif
1307 CoinFactorizationDouble pivotValue2 = region2[jRow];
1308 CoinFactorizationDouble pivotMult = pivotRegion[jRow];
1309 assert(pivotMult);
1310 CoinFactorizationDouble pivotValue2a = pivotValue2 * pivotMult;
1311 bool do2 = !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2);
1312 region2[jRow] = 0.0;
1313 CoinFactorizationDouble pivotValue1 = region1[jRow];
1314 region1[jRow] = 0.0;
1315 CoinFactorizationDouble pivotValue1a = pivotValue1 * pivotMult;
1316 bool do1 = !TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue1);
1317 if (do2) {
1318 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2a)) {
1319 region2[jRow] = pivotValue2a;
1320 index2[numberNonZeroB++] = jRow;
1321 }
1322 }
1323 if (do1) {
1324 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue1a)) {
1325 region1[jRow] = pivotValue1a;
1326 index1[numberNonZeroA++] = jRow;
1327 }
1328 }
1329 instrument_add(numberIn);
1330 if (numberIn) {
1331 if (do2) {
1332 if (!do1) {
1333 // just region 2
1334 for (CoinBigIndex j = numberIn - 1; j >= 0; j--) {
1335 CoinSimplexInt iRow = thisIndex[j];
1336 CoinFactorizationDouble value = thisElement[j];
1337 assert(value);
1338 #ifdef NO_LOAD
1339 region2[iRow] -= value * pivotValue2;
1340 #else
1341 CoinFactorizationDouble regionValue2 = region2[iRow];
1342 region2[iRow] = regionValue2 - value * pivotValue2;
1343 #endif
1344 }
1345 } else {
1346 // both
1347 instrument_add(numberIn);
1348 for (CoinBigIndex j = numberIn - 1; j >= 0; j--) {
1349 CoinSimplexInt iRow = thisIndex[j];
1350 CoinFactorizationDouble value = thisElement[j];
1351 #ifdef NO_LOAD
1352 region1[iRow] -= value * pivotValue1;
1353 region2[iRow] -= value * pivotValue2;
1354 #else
1355 CoinFactorizationDouble regionValue1 = region1[iRow];
1356 CoinFactorizationDouble regionValue2 = region2[iRow];
1357 region1[iRow] = regionValue1 - value * pivotValue1;
1358 region2[iRow] = regionValue2 - value * pivotValue2;
1359 #endif
1360 }
1361 }
1362 } else if (do1) {
1363 // just region 1
1364 for (CoinBigIndex j = numberIn - 1; j >= 0; j--) {
1365 CoinSimplexInt iRow = thisIndex[j];
1366 CoinFactorizationDouble value = thisElement[j];
1367 assert(value);
1368 #ifdef NO_LOAD
1369 region1[iRow] -= value * pivotValue1;
1370 #else
1371 CoinFactorizationDouble regionValue1 = region1[iRow];
1372 region1[iRow] = regionValue1 - value * pivotValue1;
1373 #endif
1374 }
1375 }
1376 } else {
1377 // no elements in column
1378 }
1379 }
1380 jRow = pivotLinked[jRow];
1381 }
1382 numberNonZero1 = numberNonZeroA;
1383 numberNonZero2 = numberNonZeroB;
1384 #if ABC_SMALL < 2
1385 if (factorizationStatistics()) {
1386 ftranFTCountAfterU_ += numberNonZeroA;
1387 ftranCountAfterU_ += numberNonZeroB;
1388 }
1389 #endif
1390 instrument_end();
1391 }
1392 // updateColumnU. Updates part of column (FTRANU)
1393 void CoinAbcTypeFactorization::updateColumnU(CoinIndexedVector *regionSparse
1394 #if ABC_SMALL < 2
1395 ,
1396 CoinAbcStatistics &statistics
1397 #endif
1398 #if ABC_PARALLEL
1399 ,
1400 int whichSparse
1401 #endif
1402 ) const
1403 {
1404 #if CILK_CONFLICT > 0
1405 #if ABC_PARALLEL
1406 // for conflicts
1407 #if CILK_CONFLICT > 1
1408 printf("file %s line %d which %d\n", __FILE__, __LINE__, whichSparse);
1409 #endif
1410 abc_assert((cilk_conflict & (1 << whichSparse)) == 0);
1411 cilk_conflict |= (1 << whichSparse);
1412 #else
1413 abc_assert((cilk_conflict & (1 << 0)) == 0);
1414 cilk_conflict |= (1 << 0);
1415 #endif
1416 #endif
1417 #if ABC_SMALL < 2
1418 CoinSimplexInt numberNonZero = regionSparse->getNumElements();
1419 int goSparse;
1420 // Guess at number at end
1421 if (gotSparse()) {
1422 double average = statistics.averageAfterU_ * twiddleFactor1S();
1423 assert(average);
1424 CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * average);
1425 if (newNumber < sparseThreshold_)
1426 goSparse = 1;
1427 else if (numberNonZero * 3 < numberRows_)
1428 goSparse = 0;
1429 else
1430 goSparse = -1;
1431 } else {
1432 goSparse = 0;
1433 }
1434 #endif
1435 if (!goSparse) {
1436 // densish
1437 updateColumnUDensish(regionSparse);
1438 #if ABC_SMALL < 2
1439 } else if (goSparse < 0) {
1440 // dense
1441 updateColumnUDense(regionSparse);
1442 } else {
1443 // sparse
1444 updateColumnUSparse(regionSparse
1445 #if ABC_PARALLEL
1446 ,
1447 whichSparse
1448 #endif
1449 );
1450 #endif
1451 }
1452 #if ABC_SMALL < 2
1453 if (factorizationStatistics()) {
1454 statistics.countAfterU_ += regionSparse->getNumElements();
1455 }
1456 #endif
1457 #if CILK_CONFLICT > 0
1458 #if ABC_PARALLEL
1459 // for conflicts
1460 abc_assert((cilk_conflict & (1 << whichSparse)) != 0);
1461 cilk_conflict &= ~(1 << whichSparse);
1462 #else
1463 abc_assert((cilk_conflict & (1 << 0)) != 0);
1464 cilk_conflict &= ~(1 << 0);
1465 #endif
1466 #endif
1467 }
1468 //#define COUNT_U
1469 #ifdef COUNT_U
1470 static int nUDense[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
1471 static int nUSparse[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
1472 #endif
1473 // updateColumnU. Updates part of column (FTRANU)
1474 // Updates part of column (FTRANU) real work
1475 void CoinAbcTypeFactorization::updateColumnUDensish(CoinIndexedVector *regionSparse) const
1476 {
1477 instrument_start("CoinAbcFactorizationUpdateUDensish", 2 * numberRows_);
1478 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1479 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
1480 CoinSimplexInt numberNonZero = 0;
1481 const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1482 #ifdef ABC_USE_FUNCTION_POINTERS
1483 scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1484 CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1485 #else
1486 const CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1487 const CoinBigIndex *COIN_RESTRICT startColumn = startColumnUAddress_;
1488 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1489 const CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1490 #endif
1491
1492 const CoinSimplexInt *COIN_RESTRICT pivotLinked = pivotLinkedBackwardsAddress_;
1493 CoinSimplexInt jRow = pivotLinked[numberRows_];
1494 #define ETAS_1 2
1495 #define TEST_BEFORE
1496 #ifdef TEST_BEFORE
1497 CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1498 #endif
1499 #ifdef DONT_USE_SLACKS
1500 while (jRow != -1 /*lastSlack_*/) {
1501 #else
1502 while (jRow != lastSlack_) {
1503 #endif
1504 #ifndef TEST_BEFORE
1505 CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1506 #endif
1507 if (expValue) {
1508 if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
1509 CoinFactorizationDouble pivotValue = region[jRow];
1510 #if ETAS_1 > 1
1511 CoinFactorizationDouble pivotValue2 = pivotValue * pivotRegion[jRow];
1512 #endif
1513 #ifndef ABC_USE_FUNCTION_POINTERS
1514 int number = numberInColumn[jRow];
1515 #else
1516 scatterStruct &COIN_RESTRICT scatter = scatterPointer[jRow];
1517 CoinSimplexInt number = scatter.number;
1518 #endif
1519 instrument_add(number);
1520 if (TEST_INT_NONZERO(number)) {
1521 #ifdef COUNT_U
1522 {
1523 int k = numberInColumn[jRow];
1524 if (k > 10)
1525 k = 11;
1526 nUDense[k]++;
1527 int kk = 0;
1528 for (int k = 0; k < 12; k++)
1529 kk += nUDense[k];
1530 if (kk % 1000000 == 0) {
1531 printf("ZZ");
1532 for (int k = 0; k < 12; k++)
1533 printf(" %d", nUDense[k]);
1534 printf("\n");
1535 }
1536 }
1537 #endif
1538 #ifndef ABC_USE_FUNCTION_POINTERS
1539 CoinBigIndex start = startColumn[jRow];
1540 #ifndef INLINE_IT
1541 const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
1542 const CoinSimplexInt *COIN_RESTRICT thisIndex = indexRow + start;
1543 for (CoinBigIndex j = number - 1; j >= 0; j--) {
1544 CoinSimplexInt iRow = thisIndex[j];
1545 CoinFactorizationDouble regionValue = region[iRow];
1546 CoinFactorizationDouble value = thisElement[j];
1547 assert(value);
1548 region[iRow] = regionValue - value * pivotValue;
1549 }
1550 #else
1551 CoinAbcScatterUpdate(number, pivotValue, element + start, indexRow + start, region);
1552 #endif
1553 #else
1554 CoinBigIndex start = scatter.offset;
1555 #if ABC_USE_FUNCTION_POINTERS
1556 (*(scatter.functionPointer))(number, pivotValue, elementUColumnPlus + start, region);
1557 #else
1558 CoinAbcScatterUpdate(number, pivotValue, elementUColumnPlus + start, region);
1559 #endif
1560 #endif
1561 }
1562 CoinSimplexInt kRow = jRow;
1563 jRow = pivotLinked[jRow];
1564 #ifdef TEST_BEFORE
1565 expValue = ABC_EXPONENT(region[jRow]);
1566 #endif
1567 #if ETAS_1 < 2
1568 pivotValue *= pivotRegion[kRow];
1569 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
1570 region[kRow] = pivotValue;
1571 regionIndex[numberNonZero++] = kRow;
1572 } else {
1573 region[kRow] = 0.0;
1574 }
1575 #else
1576 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2)) {
1577 region[kRow] = pivotValue2;
1578 regionIndex[numberNonZero++] = kRow;
1579 } else {
1580 region[kRow] = 0.0;
1581 }
1582 #endif
1583 } else {
1584 CoinSimplexInt kRow = jRow;
1585 jRow = pivotLinked[jRow];
1586 #ifdef TEST_BEFORE
1587 expValue = ABC_EXPONENT(region[jRow]);
1588 #endif
1589 region[kRow] = 0.0;
1590 }
1591 } else {
1592 jRow = pivotLinked[jRow];
1593 #ifdef TEST_BEFORE
1594 expValue = ABC_EXPONENT(region[jRow]);
1595 #endif
1596 }
1597 }
1598 #ifndef DONT_USE_SLACKS
1599 while (jRow != -1) {
1600 #ifndef TEST_BEFORE
1601 CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1602 #endif
1603 if (expValue) {
1604 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
1605 #if SLACK_VALUE == -1
1606 CoinFactorizationDouble pivotValue = region[jRow];
1607 assert(pivotRegion[jRow] == SLACK_VALUE);
1608 region[jRow] = -pivotValue;
1609 #endif
1610 regionIndex[numberNonZero++] = jRow;
1611 } else {
1612 region[jRow] = 0.0;
1613 }
1614 }
1615 jRow = pivotLinked[jRow];
1616 #ifdef TEST_BEFORE
1617 expValue = ABC_EXPONENT(region[jRow]);
1618 #endif
1619 }
1620 #endif
1621 regionSparse->setNumElements(numberNonZero);
1622 instrument_end();
1623 }
1624 // updateColumnU. Updates part of column (FTRANU)
1625 // Updates part of column (FTRANU) real work
1626 void CoinAbcTypeFactorization::updateColumnUDense(CoinIndexedVector *regionSparse) const
1627 {
1628 instrument_start("CoinAbcFactorizationUpdateUDensish", 2 * numberRows_);
1629 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1630 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
1631 const CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1632 CoinSimplexInt numberNonZero = 0;
1633 const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1634
1635 #ifdef ABC_USE_FUNCTION_POINTERS
1636 scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1637 CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1638 #else
1639 const CoinBigIndex *COIN_RESTRICT startColumn = startColumnUAddress_;
1640 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1641 const CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1642 #endif
1643 const CoinSimplexInt *COIN_RESTRICT pivotLinked = pivotLinkedBackwardsAddress_;
1644 CoinSimplexInt jRow = pivotLinked[numberRows_];
1645 #define ETAS_1 2
1646 #define TEST_BEFORE
1647 #ifdef TEST_BEFORE
1648 CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1649 #endif
1650 //int ixxxxxx=0;
1651 #ifdef DONT_USE_SLACKS
1652 while (jRow != -1 /*lastSlack_*/) {
1653 #else
1654 while (jRow != lastSlack_) {
1655 #endif
1656 #if 0
1657 double largest=1.0;
1658 int iLargest=-1;
1659 ixxxxxx++;
1660 for (int i=0;i<numberRows_;i++) {
1661 if (fabs(region[i])>largest) {
1662 largest=fabs(region[i]);
1663 iLargest=i;
1664 }
1665 }
1666 if (iLargest>=0)
1667 printf("largest %g on row %d after %d\n",largest,iLargest,ixxxxxx);
1668 #endif
1669 #ifndef TEST_BEFORE
1670 CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1671 #endif
1672 if (expValue) {
1673 if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
1674 CoinFactorizationDouble pivotValue = region[jRow];
1675 #if ETAS_1 > 1
1676 CoinFactorizationDouble pivotValue2 = pivotValue * pivotRegion[jRow];
1677 #endif
1678 #ifndef ABC_USE_FUNCTION_POINTERS
1679 int number = numberInColumn[jRow];
1680 #else
1681 scatterStruct &COIN_RESTRICT scatter = scatterPointer[jRow];
1682 CoinSimplexInt number = scatter.number;
1683 #endif
1684 instrument_add(number);
1685 if (TEST_INT_NONZERO(number)) {
1686 #ifdef COUNT_U
1687 {
1688 int k = numberInColumn[jRow];
1689 if (k > 10)
1690 k = 11;
1691 nUDense[k]++;
1692 int kk = 0;
1693 for (int k = 0; k < 12; k++)
1694 kk += nUDense[k];
1695 if (kk % 1000000 == 0) {
1696 printf("ZZ");
1697 for (int k = 0; k < 12; k++)
1698 printf(" %d", nUDense[k]);
1699 printf("\n");
1700 }
1701 }
1702 #endif
1703 #ifndef ABC_USE_FUNCTION_POINTERS
1704 CoinBigIndex start = startColumn[jRow];
1705 #ifndef INLINE_IT
1706 const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
1707 const CoinSimplexInt *COIN_RESTRICT thisIndex = indexRow + start;
1708 for (CoinBigIndex j = number - 1; j >= 0; j--) {
1709 CoinSimplexInt iRow = thisIndex[j];
1710 CoinFactorizationDouble regionValue = region[iRow];
1711 CoinFactorizationDouble value = thisElement[j];
1712 assert(value);
1713 region[iRow] = regionValue - value * pivotValue;
1714 }
1715 #else
1716 CoinAbcScatterUpdate(number, pivotValue, element + start, indexRow + start, region);
1717 #endif
1718 #else
1719 CoinBigIndex start = scatter.offset;
1720 #if ABC_USE_FUNCTION_POINTERS
1721 (*(scatter.functionPointer))(number, pivotValue, elementUColumnPlus + start, region);
1722 #else
1723 CoinAbcScatterUpdate(number, pivotValue, elementUColumnPlus + start, region);
1724 #endif
1725 #endif
1726 }
1727 CoinSimplexInt kRow = jRow;
1728 jRow = pivotLinked[jRow];
1729 #ifdef TEST_BEFORE
1730 expValue = ABC_EXPONENT(region[jRow]);
1731 #endif
1732 #if ETAS_1 < 2
1733 pivotValue *= pivotRegion[kRow];
1734 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
1735 region[kRow] = pivotValue;
1736 regionIndex[numberNonZero++] = kRow;
1737 } else {
1738 region[kRow] = 0.0;
1739 }
1740 #else
1741 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2)) {
1742 region[kRow] = pivotValue2;
1743 regionIndex[numberNonZero++] = kRow;
1744 } else {
1745 region[kRow] = 0.0;
1746 }
1747 #endif
1748 } else {
1749 CoinSimplexInt kRow = jRow;
1750 jRow = pivotLinked[jRow];
1751 #ifdef TEST_BEFORE
1752 expValue = ABC_EXPONENT(region[jRow]);
1753 #endif
1754 region[kRow] = 0.0;
1755 }
1756 } else {
1757 jRow = pivotLinked[jRow];
1758 #ifdef TEST_BEFORE
1759 expValue = ABC_EXPONENT(region[jRow]);
1760 #endif
1761 }
1762 }
1763 #ifndef DONT_USE_SLACKS
1764 while (jRow != -1) {
1765 #ifndef TEST_BEFORE
1766 CoinExponent expValue = ABC_EXPONENT(region[jRow]);
1767 #endif
1768 if (expValue) {
1769 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
1770 #if SLACK_VALUE == -1
1771 CoinFactorizationDouble pivotValue = region[jRow];
1772 assert(pivotRegion[jRow] == SLACK_VALUE);
1773 region[jRow] = -pivotValue;
1774 #endif
1775 regionIndex[numberNonZero++] = jRow;
1776 } else {
1777 region[jRow] = 0.0;
1778 }
1779 }
1780 jRow = pivotLinked[jRow];
1781 #ifdef TEST_BEFORE
1782 expValue = ABC_EXPONENT(region[jRow]);
1783 #endif
1784 }
1785 #endif
1786 regionSparse->setNumElements(numberNonZero);
1787 instrument_end();
1788 }
1789 #if ABC_SMALL < 2
1790 // updateColumnU. Updates part of column (FTRANU)
1791 /*
1792 Since everything is in order I should be able to do a better job of
1793 marking stuff - think. Also as L is static maybe I can do something
1794 better there (I know I could if I marked the depth of every element
1795 but that would lead to other inefficiencies.
1796 */
1797 void CoinAbcTypeFactorization::updateColumnUSparse(CoinIndexedVector *regionSparse
1798 #if ABC_PARALLEL
1799 ,
1800 int whichSparse
1801 #endif
1802 ) const
1803 {
1804 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1805 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
1806 CoinSimplexInt numberNonZero = regionSparse->getNumElements();
1807 //const CoinBigIndex * COIN_RESTRICT startColumn = startColumnUAddress_;
1808 //const CoinSimplexInt * COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1809 //const CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
1810 const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1811 // use sparse_ as temporary area
1812 // mark known to be zero
1813 #define COINFACTORIZATION_SHIFT_PER_INT2 (COINFACTORIZATION_SHIFT_PER_INT - 1)
1814 #define COINFACTORIZATION_MASK_PER_INT2 (COINFACTORIZATION_MASK_PER_INT >> 1)
1815 // need to redo if this fails (just means using CoinAbcStack to compute sizes)
1816 assert(sizeof(CoinSimplexInt) == sizeof(CoinBigIndex));
1817 CoinAbcStack *COIN_RESTRICT stackList = reinterpret_cast< CoinAbcStack * >(sparseAddress_);
1818 CoinSimplexInt *COIN_RESTRICT list = listAddress_;
1819 #ifndef ABC_USE_FUNCTION_POINTERS
1820 const CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1821 #else
1822 scatterStruct *COIN_RESTRICT scatterPointer = scatterUColumn();
1823 const CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1824 const CoinSimplexInt *COIN_RESTRICT indexRow = reinterpret_cast< const CoinSimplexInt * >(elementUColumnPlusAddress_);
1825 #endif
1826 //#define DO_CHAR2 1
1827 #if DO_CHAR2 // CHAR
1828 CoinCheckZero *COIN_RESTRICT mark = markRowAddress_;
1829 #else
1830 //BIT
1831 CoinSimplexUnsignedInt *COIN_RESTRICT mark = reinterpret_cast< CoinSimplexUnsignedInt * >(markRowAddress_);
1832 #endif
1833 #if ABC_PARALLEL
1834 //printf("PP %d %d %s\n",whichSparse,__LINE__,__FILE__);
1835 if (whichSparse) {
1836 int addAmount = whichSparse * sizeSparseArray_;
1837 stackList = reinterpret_cast< CoinAbcStack * >(reinterpret_cast< char * >(stackList) + addAmount);
1838 list = reinterpret_cast< CoinSimplexInt * >(reinterpret_cast< char * >(list) + addAmount);
1839 #if DO_CHAR2 // CHAR
1840 mark = reinterpret_cast< CoinCheckZero * >(reinterpret_cast< char * >(mark) + addAmount);
1841 #else
1842 mark = reinterpret_cast< CoinSimplexUnsignedInt * >(reinterpret_cast< char * >(mark) + addAmount);
1843 #endif
1844 }
1845 #endif
1846 #ifdef COIN_DEBUG
1847 #if DO_CHAR2 // CHAR
1848 for (CoinSimplexInt i = 0; i < maximumRowsExtra_; i++) {
1849 assert(!mark[i]);
1850 }
1851 #else
1852 //BIT
1853 for (CoinSimplexInt i = 0; i< numberRows_ >> COINFACTORIZATION_SHIFT_PER_INT; i++) {
1854 assert(!mark[i]);
1855 }
1856 #endif
1857 #endif
1858 CoinSimplexInt nList = 0;
1859 // move slacks to end of stack list
1860 int *COIN_RESTRICT putLast = list + numberRows_;
1861 int *COIN_RESTRICT put = putLast;
1862 for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
1863 CoinSimplexInt kPivot = regionIndex[i];
1864 #if DO_CHAR2 // CHAR
1865 CoinCheckZero mark_B = mark[kPivot];
1866 #else
1867 CoinSimplexUnsignedInt wordB = kPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1868 CoinSimplexUnsignedInt bitB = (kPivot & COINFACTORIZATION_MASK_PER_INT2) << 1;
1869 CoinSimplexUnsignedInt mark_B = (mark[wordB] >> bitB) & 3;
1870 #endif
1871 if (!mark_B) {
1872 #if DO_CHAR2 // CHAR
1873 mark[kPivot] = 1;
1874 #else
1875 mark[wordB] |= (1 << bitB);
1876 #endif
1877 #ifndef ABC_USE_FUNCTION_POINTERS
1878 CoinBigIndex start = startColumn[kPivot];
1879 CoinSimplexInt number = numberInColumn[kPivot];
1880 #else
1881 scatterStruct &COIN_RESTRICT scatter = scatterPointer[kPivot];
1882 CoinSimplexInt number = scatter.number;
1883 CoinBigIndex start = (scatter.offset + number) << 1;
1884 #endif
1885 stackList[0].stack = kPivot;
1886 stackList[0].next = start + number;
1887 stackList[0].start = start;
1888 CoinSimplexInt nStack = 0;
1889 while (nStack >= 0) {
1890 /* take off stack */
1891 CoinBigIndex j = stackList[nStack].next - 1;
1892 CoinBigIndex start = stackList[nStack].start;
1893 #if DO_CHAR2 == 0 // CHAR
1894 CoinSimplexUnsignedInt word0;
1895 CoinSimplexUnsignedInt bit0;
1896 #endif
1897 CoinSimplexInt jPivot;
1898 for (; j >= start; j--) {
1899 jPivot = indexRow[j];
1900 #if DO_CHAR2 // CHAR
1901 CoinCheckZero mark_j = mark[jPivot];
1902 #else
1903 word0 = jPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1904 bit0 = (jPivot & COINFACTORIZATION_MASK_PER_INT2) << 1;
1905 CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 3;
1906 #endif
1907 if (!mark_j)
1908 break;
1909 }
1910 if (j >= start) {
1911 /* put back on stack */
1912 stackList[nStack].next = j;
1913 /* and new one */
1914 #ifndef ABC_USE_FUNCTION_POINTERS
1915 CoinBigIndex start = startColumn[jPivot];
1916 CoinSimplexInt number = numberInColumn[jPivot];
1917 #else
1918 scatterStruct &COIN_RESTRICT scatter = scatterPointer[jPivot];
1919 CoinSimplexInt number = scatter.number;
1920 CoinBigIndex start = (scatter.offset + number) << 1;
1921 #endif
1922 if (number) {
1923 nStack++;
1924 stackList[nStack].stack = jPivot;
1925 stackList[nStack].next = start + number;
1926 stackList[nStack].start = start;
1927 #if DO_CHAR2 // CHAR
1928 mark[jPivot] = 1;
1929 #else
1930 mark[word0] |= (1 << bit0);
1931 #endif
1932 } else {
1933 // can do at once
1934 #ifndef NDEBUG
1935 #if DO_CHAR2 // CHAR
1936 CoinCheckZero mark_j = mark[jPivot];
1937 #else
1938 CoinSimplexUnsignedInt mark_j = (mark[word0] >> bit0) & 3;
1939 #endif
1940 assert(mark_j != 3);
1941 #endif
1942 #if ABC_SMALL < 2
1943 if (!start) {
1944 // slack - put at end
1945 --put;
1946 *put = jPivot;
1947 assert(pivotRegion[jPivot] == 1.0);
1948 } else {
1949 list[nList++] = jPivot;
1950 }
1951 #else
1952 list[nList++] = jPivot;
1953 #endif
1954 #if DO_CHAR2 // CHAR
1955 mark[jPivot] = 3;
1956 #else
1957 mark[word0] |= (3 << bit0);
1958 #endif
1959 }
1960 } else {
1961 /* finished so mark */
1962 CoinSimplexInt kPivot = stackList[nStack].stack;
1963 #if DO_CHAR2 // CHAR
1964 CoinCheckZero mark_B = mark[kPivot];
1965 #else
1966 CoinSimplexUnsignedInt wordB = kPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1967 CoinSimplexUnsignedInt bitB = (kPivot & COINFACTORIZATION_MASK_PER_INT2) << 1;
1968 CoinSimplexUnsignedInt mark_B = (mark[wordB] >> bitB) & 3;
1969 #endif
1970 assert(mark_B != 3);
1971 //if (mark_B!=3) {
1972 list[nList++] = kPivot;
1973 #if DO_CHAR2 // CHAR
1974 mark[kPivot] = 3;
1975 #else
1976 mark[wordB] |= (3 << bitB);
1977 #endif
1978 //}
1979 --nStack;
1980 }
1981 }
1982 }
1983 }
1984 instrument_start("CoinAbcFactorizationUpdateUSparse", numberNonZero + 2 * nList);
1985 numberNonZero = 0;
1986 list[-1] = -1;
1987 //assert (nList);
1988 for (CoinSimplexInt i = nList - 1; i >= 0; i--) {
1989 CoinSimplexInt iPivot = list[i];
1990 CoinExponent expValue = ABC_EXPONENT(region[iPivot]);
1991 #if DO_CHAR2 // CHAR
1992 mark[iPivot] = 0;
1993 #else
1994 CoinSimplexInt word0 = iPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
1995 mark[word0] = 0;
1996 #endif
1997 if (expValue) {
1998 if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
1999 CoinFactorizationDouble pivotValue = region[iPivot];
2000 #if ETAS_1 > 1
2001 CoinFactorizationDouble pivotValue2 = pivotValue * pivotRegion[iPivot];
2002 #endif
2003 #ifndef ABC_USE_FUNCTION_POINTERS
2004 CoinSimplexInt number = numberInColumn[iPivot];
2005 #else
2006 scatterStruct &COIN_RESTRICT scatter = scatterPointer[iPivot];
2007 CoinSimplexInt number = scatter.number;
2008 #endif
2009 if (TEST_INT_NONZERO(number)) {
2010 #ifdef COUNT_U
2011 {
2012 int k = numberInColumn[iPivot];
2013 if (k > 10)
2014 k = 11;
2015 nUSparse[k]++;
2016 int kk = 0;
2017 for (int k = 0; k < 12; k++)
2018 kk += nUSparse[k];
2019 if (kk % 1000000 == 0) {
2020 printf("ZZsparse");
2021 for (int k = 0; k < 12; k++)
2022 printf(" %d", nUSparse[k]);
2023 printf("\n");
2024 }
2025 }
2026 #endif
2027 #ifndef ABC_USE_FUNCTION_POINTERS
2028 CoinBigIndex start = startColumn[iPivot];
2029 #else
2030 //CoinBigIndex start = startColumn[iPivot];
2031 CoinBigIndex start = scatter.offset;
2032 #endif
2033 instrument_add(number);
2034 #ifndef INLINE_IT
2035 CoinBigIndex j;
2036 for (j = start; j < start + number; j++) {
2037 CoinFactorizationDouble value = element[j];
2038 assert(value);
2039 CoinSimplexInt iRow = indexRow[j];
2040 region[iRow] -= value * pivotValue;
2041 }
2042 #else
2043 #ifdef ABC_USE_FUNCTION_POINTERS
2044 #if ABC_USE_FUNCTION_POINTERS
2045 (*(scatter.functionPointer))(number, pivotValue, elementUColumnPlus + start, region);
2046 #else
2047 CoinAbcScatterUpdate(number, pivotValue, elementUColumnPlus + start, region);
2048 #endif
2049 #else
2050 CoinAbcScatterUpdate(number, pivotValue, element + start, indexRow + start, region);
2051 #endif
2052 #endif
2053 }
2054 #if ETAS_1 < 2
2055 pivotValue *= pivotRegion[iPivot];
2056 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
2057 region[iPivot] = pivotValue;
2058 regionIndex[numberNonZero++] = iPivot;
2059 }
2060 #else
2061 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue2)) {
2062 region[iPivot] = pivotValue2;
2063 regionIndex[numberNonZero++] = iPivot;
2064 } else {
2065 region[iPivot] = 0.0;
2066 }
2067 #endif
2068 } else {
2069 region[iPivot] = 0.0;
2070 }
2071 }
2072 }
2073 #if ABC_SMALL < 2
2074 // slacks
2075 for (; put < putLast; put++) {
2076 int iPivot = *put;
2077 CoinExponent expValue = ABC_EXPONENT(region[iPivot]);
2078 #if DO_CHAR2 // CHAR
2079 mark[iPivot] = 0;
2080 #else
2081 CoinSimplexInt word0 = iPivot >> COINFACTORIZATION_SHIFT_PER_INT2;
2082 mark[word0] = 0;
2083 #endif
2084 if (expValue) {
2085 if (!TEST_EXPONENT_LESS_THAN_UPDATE_TOLERANCE(expValue)) {
2086 CoinFactorizationDouble pivotValue = region[iPivot];
2087 if (!TEST_LESS_THAN_TOLERANCE_REGISTER(pivotValue)) {
2088 region[iPivot] = pivotValue;
2089 regionIndex[numberNonZero++] = iPivot;
2090 }
2091 } else {
2092 region[iPivot] = 0.0;
2093 }
2094 }
2095 }
2096 #endif
2097 #ifdef COIN_DEBUG
2098 for (CoinSimplexInt i = 0; i< maximumRowsExtra_ >> COINFACTORIZATION_SHIFT_PER_INT2; i++) {
2099 assert(!mark[i]);
2100 }
2101 #endif
2102 regionSparse->setNumElements(numberNonZero);
2103 instrument_end_and_adjust(1.3);
2104 }
2105 #endif
2106 // Store update after doing L and R - retuns false if no room
2107 bool CoinAbcTypeFactorization::storeFT(
2108 #if ABC_SMALL < 3
2109 const
2110 #endif
2111 CoinIndexedVector *updateFT)
2112 {
2113 #if ABC_SMALL < 3
2114 CoinSimplexInt numberNonZero = updateFT->getNumElements();
2115 #else
2116 CoinSimplexInt numberNonZero = numberRows_;
2117 #endif
2118 #ifndef ABC_USE_FUNCTION_POINTERS
2119 if (lengthAreaU_ >= lastEntryByColumnU_ + numberNonZero) {
2120 #else
2121 if (lengthAreaUPlus_ >= lastEntryByColumnUPlus_ + (3 * numberNonZero + 1) / 2) {
2122 scatterStruct &scatter = scatterUColumn()[numberRows_];
2123 CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
2124 #endif
2125 #if ABC_SMALL < 3
2126 const CoinFactorizationDouble *COIN_RESTRICT region = denseVector(updateFT);
2127 const CoinSimplexInt *COIN_RESTRICT regionIndex = updateFT->getIndices();
2128 #else
2129 CoinSimplexDouble *COIN_RESTRICT region = updateFT->denseVector();
2130 #endif
2131 #ifndef ABC_USE_FUNCTION_POINTERS
2132 CoinBigIndex *COIN_RESTRICT startColumnU = startColumnUAddress_;
2133 //we are going to save at end of U
2134 startColumnU[numberRows_] = lastEntryByColumnU_;
2135 CoinSimplexInt *COIN_RESTRICT putIndex = indexRowUAddress_ + lastEntryByColumnU_;
2136 CoinFactorizationDouble *COIN_RESTRICT putElement = elementUAddress_ + lastEntryByColumnU_;
2137 #else
2138 scatter.offset = lastEntryByColumnUPlus_;
2139 CoinFactorizationDouble *COIN_RESTRICT putElement = elementUColumnPlus + lastEntryByColumnUPlus_;
2140 CoinSimplexInt *COIN_RESTRICT putIndex = reinterpret_cast< CoinSimplexInt * >(putElement + numberNonZero);
2141 #endif
2142 #if ABC_SMALL < 3
2143 for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
2144 CoinSimplexInt indexValue = regionIndex[i];
2145 CoinSimplexDouble value = region[indexValue];
2146 putElement[i] = value;
2147 putIndex[i] = indexValue;
2148 }
2149 #else
2150 numberNonZero = 0;
2151 for (CoinSimplexInt i = 0; i < numberRows_; i++) {
2152 CoinExponent expValue = ABC_EXPONENT(region[i]);
2153 if (expValue) {
2154 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
2155 putElement[numberNonZero] = region[i];
2156 putIndex[numberNonZero++] = i;
2157 } else {
2158 region[i] = 0.0;
2159 }
2160 }
2161 }
2162 #endif
2163 #ifndef ABC_USE_FUNCTION_POINTERS
2164 numberInColumnAddress_[numberRows_] = numberNonZero;
2165 lastEntryByColumnU_ += numberNonZero;
2166 #else
2167 scatter.number = numberNonZero;
2168 #endif
2169 return true;
2170 } else {
2171 return false;
2172 }
2173 }
2174 CoinSimplexInt CoinAbcTypeFactorization::updateColumnFT(CoinIndexedVector ®ionSparseX)
2175 {
2176 CoinIndexedVector *regionSparse = ®ionSparseX;
2177 CoinSimplexInt numberNonZero = regionSparse->getNumElements();
2178 toLongArray(regionSparse, 0);
2179 #ifdef ABC_ORDERED_FACTORIZATION
2180 // Permute in for Ftran
2181 permuteInForFtran(regionSparseX);
2182 #endif
2183 // ******* L
2184 //printf("a\n");
2185 //regionSparse->checkClean();
2186 updateColumnL(regionSparse
2187 #if ABC_SMALL < 2
2188 ,
2189 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2190 #endif
2191 );
2192 //printf("b\n");
2193 //regionSparse->checkClean();
2194 //row bits here
2195 updateColumnR(regionSparse
2196 #if ABC_SMALL < 2
2197 ,
2198 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2199 #endif
2200 );
2201
2202 //printf("c\n");
2203 //regionSparse->checkClean();
2204 bool doFT = storeFT(regionSparse);
2205 //printf("d\n");
2206 //regionSparse->checkClean();
2207 //update counts
2208 // ******* U
2209 updateColumnU(regionSparse
2210 #if ABC_SMALL < 2
2211 ,
2212 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2213 #endif
2214 );
2215 //printf("e\n");
2216 #if ABC_DEBUG
2217 regionSparse->checkClean();
2218 #endif
2219 numberNonZero = regionSparse->getNumElements();
2220 // will be negative if no room
2221 if (doFT)
2222 return numberNonZero;
2223 else
2224 return -numberNonZero;
2225 }
2226 void CoinAbcTypeFactorization::updateColumnFT(CoinIndexedVector ®ionSparseX,
2227 CoinIndexedVector &partialUpdate,
2228 int whichSparse)
2229 {
2230 CoinIndexedVector *regionSparse = ®ionSparseX;
2231 CoinSimplexInt numberNonZero = regionSparse->getNumElements();
2232 toLongArray(regionSparse, whichSparse);
2233 #ifdef ABC_ORDERED_FACTORIZATION
2234 // Permute in for Ftran
2235 permuteInForFtran(regionSparseX);
2236 #endif
2237 // ******* L
2238 //printf("a\n");
2239 //regionSparse->checkClean();
2240 updateColumnL(regionSparse
2241 #if ABC_SMALL < 2
2242 ,
2243 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2244 #endif
2245 #if ABC_PARALLEL
2246 ,
2247 whichSparse
2248 #endif
2249 );
2250 //printf("b\n");
2251 //regionSparse->checkClean();
2252 //row bits here
2253 updateColumnR(regionSparse
2254 #if ABC_SMALL < 2
2255 ,
2256 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2257 #endif
2258 #if ABC_PARALLEL
2259 ,
2260 whichSparse
2261 #endif
2262 );
2263 numberNonZero = regionSparse->getNumElements();
2264 CoinSimplexInt *COIN_RESTRICT indexSave = partialUpdate.getIndices();
2265 CoinSimplexDouble *COIN_RESTRICT regionSave = partialUpdate.denseVector();
2266 partialUpdate.setNumElements(numberNonZero);
2267 memcpy(indexSave, regionSparse->getIndices(), numberNonZero * sizeof(CoinSimplexInt));
2268 partialUpdate.setPackedMode(false);
2269 CoinFactorizationDouble *COIN_RESTRICT region = denseVector(regionSparse);
2270 for (int i = 0; i < numberNonZero; i++) {
2271 int iRow = indexSave[i];
2272 regionSave[iRow] = region[iRow];
2273 }
2274 // ******* U
2275 updateColumnU(regionSparse
2276 #if ABC_SMALL < 2
2277 ,
2278 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2279 #endif
2280 #if ABC_PARALLEL
2281 ,
2282 whichSparse
2283 #endif
2284 );
2285 //printf("e\n");
2286 #if ABC_DEBUG
2287 regionSparse->checkClean();
2288 #endif
2289 }
2290 int CoinAbcTypeFactorization::updateColumnFTPart1(CoinIndexedVector ®ionSparse)
2291 {
2292 toLongArray(®ionSparse, 0);
2293 #ifdef ABC_ORDERED_FACTORIZATION
2294 // Permute in for Ftran
2295 permuteInForFtran(regionSparse);
2296 #endif
2297 //CoinSimplexInt numberNonZero=regionSparse.getNumElements();
2298 // ******* L
2299 //regionSparse.checkClean();
2300 updateColumnL(®ionSparse
2301 #if ABC_SMALL < 2
2302 ,
2303 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2304 #endif
2305 #if ABC_PARALLEL
2306 ,
2307 2
2308 #endif
2309 );
2310 //regionSparse.checkClean();
2311 //row bits here
2312 updateColumnR(®ionSparse
2313 #if ABC_SMALL < 2
2314 ,
2315 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2316 #endif
2317 #if ABC_PARALLEL
2318 ,
2319 2
2320 #endif
2321 );
2322
2323 //regionSparse.checkClean();
2324 bool doFT = storeFT(®ionSparse);
2325 // will be negative if no room
2326 if (doFT)
2327 return 1;
2328 else
2329 return -1;
2330 }
2331 void CoinAbcTypeFactorization::updateColumnFTPart2(CoinIndexedVector ®ionSparse)
2332 {
2333 toLongArray(®ionSparse, 0);
2334 //CoinSimplexInt numberNonZero=regionSparse.getNumElements();
2335 // ******* U
2336 updateColumnU(®ionSparse
2337 #if ABC_SMALL < 2
2338 ,
2339 reinterpret_cast< CoinAbcStatistics & >(ftranFTCountInput_)
2340 #endif
2341 #if ABC_PARALLEL
2342 ,
2343 2
2344 #endif
2345 );
2346 #if ABC_DEBUG
2347 regionSparse.checkClean();
2348 #endif
2349 }
2350 /* Updates one column (FTRAN) */
2351 void CoinAbcTypeFactorization::updateColumnCpu(CoinIndexedVector ®ionSparse, int whichCpu) const
2352 {
2353 toLongArray(®ionSparse, whichCpu);
2354 #ifdef ABC_ORDERED_FACTORIZATION
2355 // Permute in for Ftran
2356 permuteInForFtran(regionSparse);
2357 #endif
2358 // ******* L
2359 updateColumnL(®ionSparse
2360 #if ABC_SMALL < 2
2361 ,
2362 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
2363 #endif
2364 #if ABC_PARALLEL
2365 ,
2366 whichCpu
2367 #endif
2368 );
2369 //row bits here
2370 updateColumnR(®ionSparse
2371 #if ABC_SMALL < 2
2372 ,
2373 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
2374 #endif
2375 #if ABC_PARALLEL
2376 ,
2377 whichCpu
2378 #endif
2379 );
2380 //update counts
2381 // ******* U
2382 updateColumnU(®ionSparse
2383 #if ABC_SMALL < 2
2384 ,
2385 reinterpret_cast< CoinAbcStatistics & >(ftranCountInput_)
2386 #endif
2387 #if ABC_PARALLEL
2388 ,
2389 whichCpu
2390 #endif
2391 );
2392 fromLongArray(whichCpu);
2393 #if ABC_DEBUG
2394 regionSparse.checkClean();
2395 #endif
2396 }
2397 /* Updates one column (BTRAN) */
2398 void CoinAbcTypeFactorization::updateColumnTransposeCpu(CoinIndexedVector ®ionSparse, int whichCpu) const
2399 {
2400 toLongArray(®ionSparse, whichCpu);
2401
2402 CoinSimplexDouble *COIN_RESTRICT region = regionSparse.denseVector();
2403 CoinSimplexInt *COIN_RESTRICT regionIndex = regionSparse.getIndices();
2404 CoinSimplexInt numberNonZero = regionSparse.getNumElements();
2405 //#if COIN_BIG_DOUBLE==1
2406 //const CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegion_.array()+1;
2407 //#else
2408 const CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
2409 //#endif
2410
2411 #ifndef ABC_ORDERED_FACTORIZATION
2412 // can I move this
2413 #ifndef INLINE_IT3
2414 for (CoinSimplexInt i = 0; i < numberNonZero; i++) {
2415 CoinSimplexInt iRow = regionIndex[i];
2416 CoinSimplexDouble value = region[iRow];
2417 value *= pivotRegion[iRow];
2418 region[iRow] = value;
2419 }
2420 #else
2421 multiplyIndexed(numberNonZero, pivotRegion,
2422 regionIndex, region);
2423 #endif
2424 #else
2425 // Permute in for Btran
2426 permuteInForBtranAndMultiply(regionSparse);
2427 #endif
2428 // ******* U
2429 // Apply pivot region - could be combined for speed
2430 // Can only combine/move inside vector for sparse
2431 CoinSimplexInt smallestIndex = pivotLinkedForwardsAddress_[-1];
2432 #if ABC_SMALL < 2
2433 // copy of code inside transposeU
2434 bool goSparse = false;
2435 #else
2436 #define goSparse false
2437 #endif
2438 #if ABC_SMALL < 2
2439 // Guess at number at end
2440 if (gotUCopy()) {
2441 assert(btranAverageAfterU_);
2442 CoinSimplexInt newNumber = static_cast< CoinSimplexInt >(numberNonZero * btranAverageAfterU_ * twiddleBtranFactor1());
2443 if (newNumber < sparseThreshold_)
2444 goSparse = true;
2445 }
2446 #endif
2447 if (numberNonZero < 40 && (numberNonZero << 4) < numberRows_ && !goSparse) {
2448 CoinSimplexInt *COIN_RESTRICT pivotRowForward = pivotColumnAddress_;
2449 CoinSimplexInt smallest = numberRowsExtra_;
2450 for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
2451 CoinSimplexInt iRow = regionIndex[j];
2452 if (pivotRowForward[iRow] < smallest) {
2453 smallest = pivotRowForward[iRow];
2454 smallestIndex = iRow;
2455 }
2456 }
2457 }
2458 updateColumnTransposeU(®ionSparse, smallestIndex
2459 #if ABC_SMALL < 2
2460 ,
2461 reinterpret_cast< CoinAbcStatistics & >(btranCountInput_)
2462 #endif
2463 #if ABC_PARALLEL
2464 ,
2465 whichCpu
2466 #endif
2467 );
2468 //row bits here
2469 updateColumnTransposeR(®ionSparse
2470 #if ABC_SMALL < 2
2471 ,
2472 reinterpret_cast< CoinAbcStatistics & >(btranCountInput_)
2473 #endif
2474 );
2475 // ******* L
2476 updateColumnTransposeL(®ionSparse
2477 #if ABC_SMALL < 2
2478 ,
2479 reinterpret_cast< CoinAbcStatistics & >(btranCountInput_)
2480 #endif
2481 #if ABC_PARALLEL
2482 ,
2483 whichCpu
2484 #endif
2485 );
2486 #if ABC_SMALL < 3
2487 #ifdef ABC_ORDERED_FACTORIZATION
2488 // Permute out for Btran
2489 permuteOutForBtran(regionSparse);
2490 #endif
2491 #if ABC_DEBUG
2492 regionSparse.checkClean();
2493 #endif
2494 numberNonZero = regionSparse.getNumElements();
2495 #else
2496 numberNonZero = 0;
2497 for (CoinSimplexInt i = 0; i < numberRows_; i++) {
2498 CoinExponent expValue = ABC_EXPONENT(region[i]);
2499 if (expValue) {
2500 if (!TEST_EXPONENT_LESS_THAN_TOLERANCE(expValue)) {
2501 regionIndex[numberNonZero++] = i;
2502 } else {
2503 region[i] = 0.0;
2504 }
2505 }
2506 }
2507 regionSparse.setNumElements(numberNonZero);
2508 #endif
2509 }
2510 #if ABC_SMALL < 2
2511 // getRowSpaceIterate. Gets space for one Row with given length
2512 //may have to do compression (returns true)
2513 //also moves existing vector
2514 bool CoinAbcTypeFactorization::getRowSpaceIterate(CoinSimplexInt iRow,
2515 CoinSimplexInt extraNeeded)
2516 {
2517 const CoinSimplexInt *COIN_RESTRICT numberInRow = numberInRowAddress_;
2518 CoinSimplexInt number = numberInRow[iRow];
2519 CoinBigIndex *COIN_RESTRICT startRow = startRowUAddress_;
2520 CoinSimplexInt *COIN_RESTRICT indexColumnU = indexColumnUAddress_;
2521 #if CONVERTROW
2522 CoinBigIndex *COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
2523 #if CONVERTROW > 2
2524 CoinBigIndex *COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
2525 #endif
2526 #endif
2527 CoinFactorizationDouble *COIN_RESTRICT elementRowU = elementRowUAddress_;
2528 CoinBigIndex space = lengthAreaU_ - lastEntryByRowU_;
2529 CoinSimplexInt *COIN_RESTRICT nextRow = nextRowAddress_;
2530 CoinSimplexInt *COIN_RESTRICT lastRow = lastRowAddress_;
2531 if (space < extraNeeded + number + 2) {
2532 //compression
2533 CoinSimplexInt iRow = nextRow[numberRows_];
2534 CoinBigIndex put = 0;
2535 while (iRow != numberRows_) {
2536 //move
2537 CoinBigIndex get = startRow[iRow];
2538 CoinBigIndex getEnd = startRow[iRow] + numberInRow[iRow];
2539
2540 startRow[iRow] = put;
2541 CoinBigIndex i;
2542 for (i = get; i < getEnd; i++) {
2543 indexColumnU[put] = indexColumnU[i];
2544 #if CONVERTROW
2545 convertRowToColumn[put] = convertRowToColumn[i];
2546 #if CONVERTROW > 2
2547 convertColumnToRow[i] = convertColumnToRow[put];
2548 #endif
2549 #endif
2550 elementRowU[put] = elementRowU[i];
2551 put++;
2552 }
2553 iRow = nextRow[iRow];
2554 } /* endwhile */
2555 numberCompressions_++;
2556 lastEntryByRowU_ = put;
2557 space = lengthAreaU_ - put;
2558 if (space < extraNeeded + number + 2) {
2559 //need more space
2560 //if we can allocate bigger then do so and copy
2561 //if not then return so code can start again
2562 status_ = -99;
2563 return false;
2564 }
2565 }
2566 CoinBigIndex put = lastEntryByRowU_;
2567 CoinSimplexInt next = nextRow[iRow];
2568 CoinSimplexInt last = lastRow[iRow];
2569
2570 //out
2571 nextRow[last] = next;
2572 lastRow[next] = last;
2573 //in at end
2574 last = lastRow[numberRows_];
2575 nextRow[last] = iRow;
2576 lastRow[numberRows_] = iRow;
2577 lastRow[iRow] = last;
2578 nextRow[iRow] = numberRows_;
2579 //move
2580 CoinBigIndex get = startRow[iRow];
2581 startRow[iRow] = put;
2582 while (number) {
2583 number--;
2584 indexColumnU[put] = indexColumnU[get];
2585 #if CONVERTROW
2586 convertRowToColumn[put] = convertRowToColumn[get];
2587 #if CONVERTROW > 2
2588 convertColumnToRow[get] = convertColumnToRow[put];
2589 #endif
2590 #endif
2591 elementRowU[put] = elementRowU[get];
2592 put++;
2593 get++;
2594 } /* endwhile */
2595 //add four for luck
2596 lastEntryByRowU_ = put + extraNeeded + 4;
2597 return true;
2598 }
2599 #endif
2600 void CoinAbcTypeFactorization::printRegion(const CoinIndexedVector &vector, const char *where) const
2601 {
2602 //return;
2603 CoinSimplexInt numberNonZero = vector.getNumElements();
2604 //CoinSimplexInt * COIN_RESTRICT regionIndex = vector.getIndices ( );
2605 const CoinSimplexDouble *COIN_RESTRICT region = vector.denseVector();
2606 int n = 0;
2607 for (int i = 0; i < numberRows_; i++) {
2608 if (region[i])
2609 n++;
2610 }
2611 printf("==== %d nonzeros (%d in count) %s ====\n", n, numberNonZero, where);
2612 for (int i = 0; i < numberRows_; i++) {
2613 if (region[i])
2614 printf("%d %g\n", i, region[i]);
2615 }
2616 printf("==== %s ====\n", where);
2617 }
2618 void CoinAbcTypeFactorization::unpack(CoinIndexedVector *regionFrom,
2619 CoinIndexedVector *regionTo) const
2620 {
2621 // move
2622 CoinSimplexInt *COIN_RESTRICT regionIndex = regionTo->getIndices();
2623 CoinSimplexInt numberNonZero = regionFrom->getNumElements();
2624 CoinSimplexInt *COIN_RESTRICT index = regionFrom->getIndices();
2625 CoinSimplexDouble *COIN_RESTRICT region = regionTo->denseVector();
2626 CoinSimplexDouble *COIN_RESTRICT array = regionFrom->denseVector();
2627 for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
2628 CoinSimplexInt iRow = index[j];
2629 CoinSimplexDouble value = array[j];
2630 array[j] = 0.0;
2631 region[iRow] = value;
2632 regionIndex[j] = iRow;
2633 }
2634 regionTo->setNumElements(numberNonZero);
2635 }
2636 void CoinAbcTypeFactorization::pack(CoinIndexedVector *regionFrom,
2637 CoinIndexedVector *regionTo) const
2638 {
2639 // move
2640 CoinSimplexInt *COIN_RESTRICT regionIndex = regionFrom->getIndices();
2641 CoinSimplexInt numberNonZero = regionFrom->getNumElements();
2642 CoinSimplexInt *COIN_RESTRICT index = regionTo->getIndices();
2643 CoinSimplexDouble *COIN_RESTRICT region = regionFrom->denseVector();
2644 CoinSimplexDouble *COIN_RESTRICT array = regionTo->denseVector();
2645 for (CoinSimplexInt j = 0; j < numberNonZero; j++) {
2646 CoinSimplexInt iRow = regionIndex[j];
2647 CoinSimplexDouble value = region[iRow];
2648 region[iRow] = 0.0;
2649 array[j] = value;
2650 index[j] = iRow;
2651 }
2652 regionTo->setNumElements(numberNonZero);
2653 }
2654 // Set up addresses from arrays
2655 void CoinAbcTypeFactorization::doAddresses()
2656 {
2657 pivotColumnAddress_ = pivotColumn_.array();
2658 permuteAddress_ = permute_.array();
2659 pivotRegionAddress_ = pivotRegion_.array() + 1;
2660 elementUAddress_ = elementU_.array();
2661 indexRowUAddress_ = indexRowU_.array();
2662 numberInColumnAddress_ = numberInColumn_.array();
2663 numberInColumnPlusAddress_ = numberInColumnPlus_.array();
2664 #ifdef ABC_USE_FUNCTION_POINTERS
2665 scatterPointersUColumnAddress_ = reinterpret_cast< scatterStruct * >(scatterUColumn_.array());
2666 #endif
2667 startColumnUAddress_ = startColumnU_.array();
2668 #if CONVERTROW
2669 convertRowToColumnUAddress_ = convertRowToColumnU_.array();
2670 #if CONVERTROW > 1
2671 convertColumnToRowUAddress_ = convertColumnToRowU_.array();
2672 #endif
2673 #endif
2674 #if ABC_SMALL < 2
2675 elementRowUAddress_ = elementRowU_.array();
2676 #endif
2677 startRowUAddress_ = startRowU_.array();
2678 numberInRowAddress_ = numberInRow_.array();
2679 indexColumnUAddress_ = indexColumnU_.array();
2680 firstCountAddress_ = firstCount_.array();
2681 nextCountAddress_ = nextCount();
2682 lastCountAddress_ = lastCount();
2683 nextColumnAddress_ = nextColumn_.array();
2684 lastColumnAddress_ = lastColumn_.array();
2685 nextRowAddress_ = nextRow_.array();
2686 lastRowAddress_ = lastRow_.array();
2687 saveColumnAddress_ = saveColumn_.array();
2688 //saveColumnAddress2_ = saveColumn_.array()+numberRows_;
2689 elementLAddress_ = elementL_.array();
2690 indexRowLAddress_ = indexRowL_.array();
2691 startColumnLAddress_ = startColumnL_.array();
2692 #if ABC_SMALL < 2
2693 startRowLAddress_ = startRowL_.array();
2694 #endif
2695 pivotLinkedBackwardsAddress_ = pivotLinkedBackwards();
2696 pivotLinkedForwardsAddress_ = pivotLinkedForwards();
2697 pivotLOrderAddress_ = pivotLOrder();
2698 startColumnRAddress_ = startColumnR();
2699 // next two are recomputed cleanup
2700 elementRAddress_ = elementL_.array() + lengthL_;
2701 indexRowRAddress_ = indexRowL_.array() + lengthL_;
2702 #if ABC_SMALL < 2
2703 indexColumnLAddress_ = indexColumnL_.array();
2704 elementByRowLAddress_ = elementByRowL_.array();
2705 #endif
2706 #if ABC_DENSE_CODE
2707 denseAreaAddress_ = denseArea_.array();
2708 #endif
2709 workAreaAddress_ = workArea_.array();
2710 workArea2Address_ = workArea2_.array();
2711 #if ABC_SMALL < 2
2712 sparseAddress_ = sparse_.array();
2713 CoinAbcStack *stackList = reinterpret_cast< CoinAbcStack * >(sparseAddress_);
2714 listAddress_ = reinterpret_cast< CoinSimplexInt * >(stackList + numberRows_);
2715 markRowAddress_ = reinterpret_cast< CoinCheckZero * >(listAddress_ + numberRows_);
2716 #endif
2717 }
2718 #endif
2719
2720 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2721 */
2722