1 /* $Id: CoinFactorization3.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #if defined(_MSC_VER)
7 // Turn off compiler warning about long names
8 #pragma warning(disable : 4786)
9 #endif
10
11 #include "CoinUtilsConfig.h"
12
13 #include <cassert>
14 #include <cstdio>
15
16 #include "CoinFactorization.hpp"
17 #include "CoinIndexedVector.hpp"
18 #include "CoinHelperFunctions.hpp"
19 #include "CoinTime.hpp"
20 #include <stdio.h>
21 #include <iostream>
22 #if COIN_FACTORIZATION_DENSE_CODE == 1
23 // using simple lapack interface
24 extern "C" {
25 /** LAPACK Fortran subroutine DGETRS. */
26 void F77_FUNC(dgetrs, DGETRS)(char *trans, cipfint *n,
27 cipfint *nrhs, const double *A, cipfint *ldA,
28 cipfint *ipiv, double *B, cipfint *ldB, ipfint *info,
29 int trans_len);
30 }
31 #elif COIN_FACTORIZATION_DENSE_CODE == 2
32 // C interface
33 enum CBLAS_ORDER { CblasRowMajor = 101,
34 CblasColMajor = 102 };
35 enum CBLAS_TRANSPOSE { CblasNoTrans = 111,
36 CblasTrans = 112 };
37 extern "C" {
38 int clapack_dgetrs(const enum CBLAS_ORDER Order,
39 const enum CBLAS_TRANSPOSE Trans,
40 const int N, const int NRHS,
41 const double *A, const int lda, const int *ipiv, double *B,
42 const int ldb);
43 }
44 #elif COIN_FACTORIZATION_DENSE_CODE == 3
45 // Intel compiler
46 #include "mkl_lapacke.h"
47 #endif
48 // For semi-sparse
49 #define BITS_PER_CHECK 8
50 #define CHECK_SHIFT 3
51 typedef unsigned char CoinCheckZero;
52 #ifdef CLP_FACTORIZATION_INSTRUMENT
53 extern double externalTimeStart;
54 extern double timeInFactorize;
55 extern double timeInUpdate;
56 extern double timeInFactorizeFake;
57 extern double timeInUpdateFake1;
58 extern double timeInUpdateFake2;
59 extern double timeInUpdateTranspose;
60 extern double timeInUpdateFT;
61 extern double timeInUpdateTwoFT;
62 extern double timeInReplace;
63 extern int numberUpdate;
64 extern int numberUpdateTranspose;
65 extern int numberUpdateFT;
66 extern int numberUpdateTwoFT;
67 extern int numberReplace;
68 extern int currentLengthR;
69 extern int currentLengthU;
70 extern int currentTakeoutU;
71 extern double averageLengthR;
72 extern double averageLengthL;
73 extern double averageLengthU;
74 extern double scaledLengthDense;
75 extern double scaledLengthDenseSquared;
76 extern double scaledLengthL;
77 extern double scaledLengthR;
78 extern double scaledLengthU;
79 extern int startLengthU;
80 extern int endLengthU;
81 #endif
82
83 //:class CoinFactorization. Deals with Factorization and Updates
84
85 /* Updates one column (FTRAN) from region2 and permutes.
86 region1 starts as zero
87 Note - if regionSparse2 packed on input - will be packed on output
88 - returns un-permuted result in region2 and region1 is zero */
updateColumn(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2,bool noPermute) const89 int CoinFactorization::updateColumn(CoinIndexedVector *regionSparse,
90 CoinIndexedVector *regionSparse2,
91 bool noPermute)
92 const
93 {
94 #ifdef CLP_FACTORIZATION_INSTRUMENT
95 double startTimeX = CoinCpuTime();
96 #endif
97 //permute and move indices into index array
98 int *COIN_RESTRICT regionIndex = regionSparse->getIndices();
99 int numberNonZero;
100 const int *permute = permute_.array();
101 double *COIN_RESTRICT region = regionSparse->denseVector();
102
103 #ifndef CLP_FACTORIZATION
104 if (!noPermute) {
105 #endif
106 numberNonZero = regionSparse2->getNumElements();
107 int *COIN_RESTRICT index = regionSparse2->getIndices();
108 double *COIN_RESTRICT array = regionSparse2->denseVector();
109 #ifndef CLP_FACTORIZATION
110 bool packed = regionSparse2->packedMode();
111 if (packed) {
112 for (int j = 0; j < numberNonZero; j++) {
113 int iRow = index[j];
114 double value = array[j];
115 array[j] = 0.0;
116 iRow = permute[iRow];
117 region[iRow] = value;
118 regionIndex[j] = iRow;
119 }
120 } else {
121 #else
122 assert(!regionSparse2->packedMode());
123 #endif
124 for (int j = 0; j < numberNonZero; j++) {
125 int iRow = index[j];
126 double value = array[iRow];
127 array[iRow] = 0.0;
128 iRow = permute[iRow];
129 region[iRow] = value;
130 regionIndex[j] = iRow;
131 }
132 #ifndef CLP_FACTORIZATION
133 }
134 #endif
135 regionSparse->setNumElements(numberNonZero);
136 #ifndef CLP_FACTORIZATION
137 } else {
138 numberNonZero = regionSparse->getNumElements();
139 }
140 #endif
141 if (collectStatistics_) {
142 numberFtranCounts_++;
143 ftranCountInput_ += numberNonZero;
144 }
145
146 // ******* L
147 updateColumnL(regionSparse, regionIndex);
148 if (collectStatistics_)
149 ftranCountAfterL_ += regionSparse->getNumElements();
150 //permute extra
151 //row bits here
152 updateColumnR(regionSparse);
153 if (collectStatistics_)
154 ftranCountAfterR_ += regionSparse->getNumElements();
155
156 //update counts
157 // ******* U
158 updateColumnU(regionSparse, regionIndex);
159 if (!doForrestTomlin_) {
160 // Do PFI after everything else
161 updateColumnPFI(regionSparse);
162 }
163 #ifdef CLP_FACTORIZATION_INSTRUMENT
164 numberUpdate++;
165 timeInUpdate += CoinCpuTime() - startTimeX;
166 averageLengthR += lengthR_;
167 averageLengthU += lengthU_;
168 averageLengthL += lengthL_;
169 #endif
170 if (!noPermute) {
171 permuteBack(regionSparse, regionSparse2);
172 return regionSparse2->getNumElements();
173 } else {
174 return regionSparse->getNumElements();
175 }
176 }
177 // Permutes back at end of updateColumn
permuteBack(CoinIndexedVector * regionSparse,CoinIndexedVector * outVector) const178 void CoinFactorization::permuteBack(CoinIndexedVector *regionSparse,
179 CoinIndexedVector *outVector) const
180 {
181 // permute back
182 int oldNumber = regionSparse->getNumElements();
183 const int *COIN_RESTRICT regionIndex = regionSparse->getIndices();
184 double *COIN_RESTRICT region = regionSparse->denseVector();
185 int *COIN_RESTRICT outIndex = outVector->getIndices();
186 double *COIN_RESTRICT out = outVector->denseVector();
187 const int *COIN_RESTRICT permuteBack = pivotColumnBack();
188 int number = 0;
189 if (outVector->packedMode()) {
190 for (int j = 0; j < oldNumber; j++) {
191 int iRow = regionIndex[j];
192 double value = region[iRow];
193 region[iRow] = 0.0;
194 if (fabs(value) > zeroTolerance_) {
195 iRow = permuteBack[iRow];
196 outIndex[number] = iRow;
197 out[number++] = value;
198 }
199 }
200 } else {
201 int j = 0;
202 if ((oldNumber & 1) != 0) {
203 int iRow = regionIndex[0];
204 j++;
205 double value = region[iRow];
206 region[iRow] = 0.0;
207 if (fabs(value) > zeroTolerance_) {
208 iRow = permuteBack[iRow];
209 outIndex[number++] = iRow;
210 out[iRow] = value;
211 }
212 }
213 for (; j < oldNumber; j += 2) {
214 int iRow0 = regionIndex[j];
215 int iRow1 = regionIndex[j + 1];
216 double value0 = region[iRow0];
217 bool good0 = fabs(value0) > zeroTolerance_;
218 double value1 = region[iRow1];
219 bool good1 = fabs(value1) > zeroTolerance_;
220 region[iRow0] = 0.0;
221 region[iRow1] = 0.0;
222 if (good0) {
223 iRow0 = permuteBack[iRow0];
224 outIndex[number++] = iRow0;
225 out[iRow0] = value0;
226 }
227 if (good1) {
228 iRow1 = permuteBack[iRow1];
229 outIndex[number++] = iRow1;
230 out[iRow1] = value1;
231 }
232 }
233 }
234 outVector->setNumElements(number);
235 regionSparse->setNumElements(0);
236 }
237 // updateColumnL. Updates part of column (FTRANL)
updateColumnL(CoinIndexedVector * regionSparse,int * COIN_RESTRICT regionIndex) const238 void CoinFactorization::updateColumnL(CoinIndexedVector *regionSparse,
239 int *COIN_RESTRICT regionIndex) const
240 {
241 if (numberL_) {
242 int number = regionSparse->getNumElements();
243 int goSparse;
244 // Guess at number at end
245 if (sparseThreshold_ > 0) {
246 if (ftranAverageAfterL_) {
247 int newNumber = static_cast< int >(number * ftranAverageAfterL_);
248 if (newNumber < sparseThreshold_ && (numberL_ << 2) > newNumber)
249 goSparse = 2;
250 else if (newNumber < sparseThreshold2_ && (numberL_ << 1) > newNumber)
251 goSparse = 1;
252 else
253 goSparse = 0;
254 } else {
255 if (number < sparseThreshold_ && (numberL_ << 2) > number)
256 goSparse = 2;
257 else
258 goSparse = 0;
259 }
260 } else {
261 goSparse = 0;
262 }
263 switch (goSparse) {
264 case 0: // densish
265 updateColumnLDensish(regionSparse, regionIndex);
266 break;
267 case 1: // middling
268 updateColumnLSparsish(regionSparse, regionIndex);
269 break;
270 case 2: // sparse
271 updateColumnLSparse(regionSparse, regionIndex);
272 break;
273 }
274 }
275 #ifdef COIN_FACTORIZATION_DENSE_CODE
276 if (numberDense_) {
277 //take off list
278 int lastSparse = numberRows_ - numberDense_;
279 int number = regionSparse->getNumElements();
280 double *COIN_RESTRICT region = regionSparse->denseVector();
281 int i = 0;
282 bool doDense = false;
283 while (i < number) {
284 int iRow = regionIndex[i];
285 if (iRow >= lastSparse) {
286 doDense = true;
287 regionIndex[i] = regionIndex[--number];
288 } else {
289 i++;
290 }
291 }
292 if (doDense) {
293 #if COIN_FACTORIZATION_DENSE_CODE == 1
294 char trans = 'N';
295 int ione = 1;
296 int info;
297 F77_FUNC(dgetrs, DGETRS)
298 (&trans, &numberDense_, &ione, denseAreaAddress_, &numberDense_,
299 densePermute_, region + lastSparse, &numberDense_, &info, 1);
300 #elif COIN_FACTORIZATION_DENSE_CODE == 2
301 clapack_dgetrs(CblasColMajor, CblasNoTrans, numberDense_, 1,
302 denseAreaAddress_, numberDense_, densePermute_,
303 region + lastSparse, numberDense_);
304 #elif COIN_FACTORIZATION_DENSE_CODE == 3
305 LAPACKE_dgetrs(LAPACK_COL_MAJOR, 'N', numberDense_, 1,
306 denseAreaAddress_, numberDense_, densePermute_,
307 region + lastSparse, numberDense_);
308 #endif
309 for (int i = lastSparse; i < numberRows_; i++) {
310 double value = region[i];
311 if (value) {
312 if (fabs(value) >= 1.0e-15)
313 regionIndex[number++] = i;
314 else
315 region[i] = 0.0;
316 }
317 }
318 regionSparse->setNumElements(number);
319 }
320 }
321 #endif
322 }
323 // Updates part of column (FTRANL) when densish
updateColumnLDensish(CoinIndexedVector * regionSparse,int * COIN_RESTRICT regionIndex) const324 void CoinFactorization::updateColumnLDensish(CoinIndexedVector *regionSparse,
325 int *COIN_RESTRICT regionIndex)
326 const
327 {
328 double *COIN_RESTRICT region = regionSparse->denseVector();
329 int number = regionSparse->getNumElements();
330 int numberNonZero;
331 double tolerance = zeroTolerance_;
332
333 numberNonZero = 0;
334
335 const int *COIN_RESTRICT startColumn = startColumnL_.array();
336 const int *COIN_RESTRICT indexRow = indexRowL_.array();
337 const CoinFactorizationDouble *COIN_RESTRICT element = elementL_.array();
338 int last = numberRows_;
339 assert(last == baseL_ + numberL_);
340 #if COIN_FACTORIZATION_DENSE_CODE
341 //can take out last bit of sparse L as empty
342 last -= numberDense_;
343 #endif
344 int smallestIndex = numberRowsExtra_;
345 // do easy ones
346 for (int k = 0; k < number; k++) {
347 int iPivot = regionIndex[k];
348 if (iPivot >= baseL_)
349 smallestIndex = CoinMin(iPivot, smallestIndex);
350 else
351 regionIndex[numberNonZero++] = iPivot;
352 }
353 // now others
354 for (int i = smallestIndex; i < last; i++) {
355 CoinFactorizationDouble pivotValue = region[i];
356
357 if (fabs(pivotValue) > tolerance) {
358 int start = startColumn[i];
359 int end = startColumn[i + 1];
360 for (int j = start; j < end; j++) {
361 int iRow = indexRow[j];
362 CoinFactorizationDouble result = region[iRow];
363 CoinFactorizationDouble value = element[j];
364
365 region[iRow] = result - value * pivotValue;
366 }
367 regionIndex[numberNonZero++] = i;
368 } else {
369 region[i] = 0.0;
370 }
371 }
372 // and dense
373 for (int i = last; i < numberRows_; i++) {
374 CoinFactorizationDouble pivotValue = region[i];
375 if (fabs(pivotValue) > tolerance) {
376 regionIndex[numberNonZero++] = i;
377 } else {
378 region[i] = 0.0;
379 }
380 }
381 regionSparse->setNumElements(numberNonZero);
382 }
383 // Updates part of column (FTRANL) when sparsish
updateColumnLSparsish(CoinIndexedVector * regionSparse,int * COIN_RESTRICT regionIndex) const384 void CoinFactorization::updateColumnLSparsish(CoinIndexedVector *regionSparse,
385 int *COIN_RESTRICT regionIndex)
386 const
387 {
388 double *COIN_RESTRICT region = regionSparse->denseVector();
389 int number = regionSparse->getNumElements();
390 int numberNonZero;
391 double tolerance = zeroTolerance_;
392
393 numberNonZero = 0;
394
395 const int *startColumn = startColumnL_.array();
396 const int *indexRow = indexRowL_.array();
397 const CoinFactorizationDouble *element = elementL_.array();
398 int last = numberRows_;
399 assert(last == baseL_ + numberL_);
400 #if COIN_FACTORIZATION_DENSE_CODE
401 //can take out last bit of sparse L as empty
402 last -= numberDense_;
403 #endif
404 // mark known to be zero
405 int nInBig = sizeof(int) / sizeof(int);
406 #if ABOCA_LITE_FACTORIZATION == 0
407 #define sparseOffset 0
408 #else
409 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
410 assert(!sparseOffset);
411 #endif
412 CoinCheckZero *COIN_RESTRICT mark = reinterpret_cast< CoinCheckZero * >(sparse_.array() + (2 + nInBig) * maximumRowsExtra_ + sparseOffset);
413 int smallestIndex = numberRowsExtra_;
414 // do easy ones
415 for (int k = 0; k < number; k++) {
416 int iPivot = regionIndex[k];
417 if (iPivot < baseL_) {
418 regionIndex[numberNonZero++] = iPivot;
419 } else {
420 smallestIndex = CoinMin(iPivot, smallestIndex);
421 int iWord = iPivot >> CHECK_SHIFT;
422 int iBit = iPivot - (iWord << CHECK_SHIFT);
423 if (mark[iWord]) {
424 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
425 } else {
426 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
427 }
428 }
429 }
430 // now others
431 // First do up to convenient power of 2
432 int jLast = (smallestIndex + BITS_PER_CHECK - 1) >> CHECK_SHIFT;
433 jLast = CoinMin((jLast << CHECK_SHIFT), last);
434 int i;
435 for (i = smallestIndex; i < jLast; i++) {
436 CoinFactorizationDouble pivotValue = region[i];
437 int start = startColumn[i];
438 int end = startColumn[i + 1];
439
440 if (fabs(pivotValue) > tolerance) {
441 for (int j = start; j < end; j++) {
442 int iRow = indexRow[j];
443 CoinFactorizationDouble result = region[iRow];
444 CoinFactorizationDouble value = element[j];
445 region[iRow] = result - value * pivotValue;
446 int iWord = iRow >> CHECK_SHIFT;
447 int iBit = iRow - (iWord << CHECK_SHIFT);
448 if (mark[iWord]) {
449 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
450 } else {
451 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
452 }
453 }
454 regionIndex[numberNonZero++] = i;
455 } else {
456 region[i] = 0.0;
457 }
458 }
459
460 int kLast = last >> CHECK_SHIFT;
461 if (jLast < last) {
462 // now do in chunks
463 for (int k = (jLast >> CHECK_SHIFT); k < kLast; k++) {
464 unsigned int iMark = mark[k];
465 if (iMark) {
466 // something in chunk - do all (as imark may change)
467 i = k << CHECK_SHIFT;
468 int iLast = i + BITS_PER_CHECK;
469 for (; i < iLast; i++) {
470 CoinFactorizationDouble pivotValue = region[i];
471 int start = startColumn[i];
472 int end = startColumn[i + 1];
473
474 if (fabs(pivotValue) > tolerance) {
475 int j;
476 for (j = start; j < end; j++) {
477 int iRow = indexRow[j];
478 CoinFactorizationDouble result = region[iRow];
479 CoinFactorizationDouble value = element[j];
480 region[iRow] = result - value * pivotValue;
481 int iWord = iRow >> CHECK_SHIFT;
482 int iBit = iRow - (iWord << CHECK_SHIFT);
483 if (mark[iWord]) {
484 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
485 } else {
486 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
487 }
488 }
489 regionIndex[numberNonZero++] = i;
490 } else {
491 region[i] = 0.0;
492 }
493 }
494 mark[k] = 0; // zero out marked
495 }
496 }
497 i = kLast << CHECK_SHIFT;
498 }
499 for (; i < last; i++) {
500 CoinFactorizationDouble pivotValue = region[i];
501 int start = startColumn[i];
502 int end = startColumn[i + 1];
503
504 if (fabs(pivotValue) > tolerance) {
505 for (int j = start; j < end; j++) {
506 int iRow = indexRow[j];
507 CoinFactorizationDouble result = region[iRow];
508 CoinFactorizationDouble value = element[j];
509 region[iRow] = result - value * pivotValue;
510 }
511 regionIndex[numberNonZero++] = i;
512 } else {
513 region[i] = 0.0;
514 }
515 }
516 // Now dense part
517 for (; i < numberRows_; i++) {
518 double pivotValue = region[i];
519 if (fabs(pivotValue) > tolerance) {
520 regionIndex[numberNonZero++] = i;
521 } else {
522 region[i] = 0.0;
523 }
524 }
525 // zero out ones that might have been skipped
526 mark[smallestIndex >> CHECK_SHIFT] = 0;
527 int kkLast = (numberRows_ + BITS_PER_CHECK - 1) >> CHECK_SHIFT;
528 CoinZeroN(mark + kLast, kkLast - kLast);
529 regionSparse->setNumElements(numberNonZero);
530 }
531 // Updates part of column (FTRANL) when sparse
updateColumnLSparse(CoinIndexedVector * regionSparse,int * COIN_RESTRICT regionIndex) const532 void CoinFactorization::updateColumnLSparse(CoinIndexedVector *regionSparse,
533 int *COIN_RESTRICT regionIndex)
534 const
535 {
536 double *COIN_RESTRICT region = regionSparse->denseVector();
537 int number = regionSparse->getNumElements();
538 int numberNonZero;
539 double tolerance = zeroTolerance_;
540
541 numberNonZero = 0;
542
543 const int *startColumn = startColumnL_.array();
544 const int *indexRow = indexRowL_.array();
545 const CoinFactorizationDouble *element = elementL_.array();
546 // use sparse_ as temporary area
547 // mark known to be zero
548 #if ABOCA_LITE_FACTORIZATION == 0
549 #define sparseOffset 0
550 #else
551 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
552 assert(!sparseOffset);
553 #endif
554 int *COIN_RESTRICT stack = sparse_.array() + sparseOffset; /* pivot */
555 int *COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
556 int *COIN_RESTRICT next = reinterpret_cast< int * >(list + maximumRowsExtra_); /* jnext */
557 char *COIN_RESTRICT mark = reinterpret_cast< char * >(next + maximumRowsExtra_);
558 int nList;
559 #ifdef COIN_DEBUG
560 for (int i = 0; i < maximumRowsExtra_; i++) {
561 assert(!mark[i]);
562 }
563 #endif
564 nList = 0;
565 for (int k = 0; k < number; k++) {
566 int kPivot = regionIndex[k];
567 if (kPivot >= baseL_) {
568 assert(kPivot < numberRowsExtra_);
569 //if (kPivot>=numberRowsExtra_) abort();
570 if (!mark[kPivot]) {
571 stack[0] = kPivot;
572 int j = startColumn[kPivot + 1] - 1;
573 int nStack = 0;
574 while (nStack >= 0) {
575 /* take off stack */
576 if (j >= startColumn[kPivot]) {
577 int jPivot = indexRow[j--];
578 assert(jPivot >= baseL_ && jPivot < numberRowsExtra_);
579 //if (jPivot<baseL_||jPivot>=numberRowsExtra_) abort();
580 /* put back on stack */
581 next[nStack] = j;
582 if (!mark[jPivot]) {
583 /* and new one */
584 kPivot = jPivot;
585 j = startColumn[kPivot + 1] - 1;
586 stack[++nStack] = kPivot;
587 assert(kPivot < numberRowsExtra_);
588 //if (kPivot>=numberRowsExtra_) abort();
589 mark[kPivot] = 1;
590 next[nStack] = j;
591 }
592 } else {
593 /* finished so mark */
594 list[nList++] = kPivot;
595 mark[kPivot] = 1;
596 --nStack;
597 if (nStack >= 0) {
598 kPivot = stack[nStack];
599 assert(kPivot < numberRowsExtra_);
600 j = next[nStack];
601 }
602 }
603 }
604 }
605 } else {
606 // just put on list
607 regionIndex[numberNonZero++] = kPivot;
608 }
609 }
610 for (int i = nList - 1; i >= 0; i--) {
611 int iPivot = list[i];
612 mark[iPivot] = 0;
613 CoinFactorizationDouble pivotValue = region[iPivot];
614 if (fabs(pivotValue) > tolerance) {
615 regionIndex[numberNonZero++] = iPivot;
616 for (int j = startColumn[iPivot];
617 j < startColumn[iPivot + 1]; j++) {
618 int iRow = indexRow[j];
619 CoinFactorizationDouble value = element[j];
620 region[iRow] -= value * pivotValue;
621 }
622 } else {
623 region[iPivot] = 0.0;
624 }
625 }
626 regionSparse->setNumElements(numberNonZero);
627 }
628 /* Updates one column (FTRAN) from region2
629 Tries to do FT update
630 number returned is negative if no room.
631 Also updates region3
632 region1 starts as zero and is zero at end */
updateTwoColumnsFT(CoinIndexedVector * regionSparse1,CoinIndexedVector * regionSparse2,CoinIndexedVector * regionSparse3,bool noPermuteRegion3)633 int CoinFactorization::updateTwoColumnsFT(CoinIndexedVector *regionSparse1,
634 CoinIndexedVector *regionSparse2,
635 CoinIndexedVector *regionSparse3,
636 bool noPermuteRegion3)
637 {
638 #ifdef CLP_FACTORIZATION_INSTRUMENT
639 double startTimeX = CoinCpuTime();
640 #endif
641 #if 1
642 //#ifdef NDEBUG
643 //#undef NDEBUG
644 //#endif
645 //#define COIN_DEBUG
646 #ifdef COIN_DEBUG
647 regionSparse1->checkClean();
648 CoinIndexedVector save2(*regionSparse2);
649 CoinIndexedVector save3(*regionSparse3);
650 #endif
651 CoinIndexedVector *regionFT;
652 CoinIndexedVector *regionUpdate;
653 int *COIN_RESTRICT regionIndex;
654 int numberNonZero;
655 const int *permute = permute_.array();
656 int *COIN_RESTRICT index;
657 double *COIN_RESTRICT region;
658 if (!noPermuteRegion3) {
659 regionFT = regionSparse3;
660 regionUpdate = regionSparse1;
661 //permute and move indices into index array
662 regionIndex = regionUpdate->getIndices();
663 //int numberNonZero;
664 region = regionUpdate->denseVector();
665
666 numberNonZero = regionSparse3->getNumElements();
667 int *COIN_RESTRICT index = regionSparse3->getIndices();
668 double *COIN_RESTRICT array = regionSparse3->denseVector();
669 assert(!regionSparse3->packedMode());
670 for (int j = 0; j < numberNonZero; j++) {
671 int iRow = index[j];
672 double value = array[iRow];
673 array[iRow] = 0.0;
674 iRow = permute[iRow];
675 region[iRow] = value;
676 regionIndex[j] = iRow;
677 }
678 regionUpdate->setNumElements(numberNonZero);
679 } else {
680 regionFT = regionSparse1;
681 regionUpdate = regionSparse3;
682 }
683 //permute and move indices into index array (in U)
684 regionIndex = regionFT->getIndices();
685 numberNonZero = regionSparse2->getNumElements();
686 index = regionSparse2->getIndices();
687 region = regionFT->denseVector();
688 double *COIN_RESTRICT array = regionSparse2->denseVector();
689 int *COIN_RESTRICT startColumnU = startColumnU_.array();
690 int start = startColumnU[maximumColumnsExtra_];
691 startColumnU[numberColumnsExtra_] = start;
692 regionIndex = indexRowU_.array() + start;
693
694 #ifndef ABC_USE_COIN_FACTORIZATION
695 assert(regionSparse2->packedMode());
696 #else
697 if (regionSparse2->packedMode()) {
698 #endif
699 for (int j = 0; j < numberNonZero; j++) {
700 int iRow = index[j];
701 double value = array[j];
702 array[j] = 0.0;
703 iRow = permute[iRow];
704 region[iRow] = value;
705 regionIndex[j] = iRow;
706 }
707 #ifdef ABC_USE_COIN_FACTORIZATION
708 }
709 else
710 {
711 // not packed
712 for (int j = 0; j < numberNonZero; j++) {
713 int iRow = index[j];
714 double value = array[iRow];
715 array[iRow] = 0.0;
716 iRow = permute[iRow];
717 region[iRow] = value;
718 regionIndex[j] = iRow;
719 }
720 }
721 #endif
722 regionFT->setNumElements(numberNonZero);
723 if (collectStatistics_) {
724 numberFtranCounts_ += 2;
725 ftranCountInput_ += regionFT->getNumElements() + regionUpdate->getNumElements();
726 }
727
728 // ******* L
729 updateColumnL(regionFT, regionIndex);
730 updateColumnL(regionUpdate, regionUpdate->getIndices());
731 if (collectStatistics_)
732 ftranCountAfterL_ += regionFT->getNumElements() + regionUpdate->getNumElements();
733 //permute extra
734 //row bits here
735 updateColumnRFT(regionFT, regionIndex);
736 updateColumnR(regionUpdate);
737 if (collectStatistics_)
738 ftranCountAfterR_ += regionFT->getNumElements() + regionUpdate->getNumElements();
739 // ******* U - see if densish
740 // Guess at number at end
741 int goSparse = 0;
742 if (sparseThreshold_ > 0) {
743 int numberNonZero = (regionUpdate->getNumElements() + regionFT->getNumElements()) >> 1;
744 if (ftranAverageAfterR_) {
745 int newNumber = static_cast< int >(numberNonZero * ftranAverageAfterU_);
746 if (newNumber < sparseThreshold_)
747 goSparse = 2;
748 else if (newNumber < sparseThreshold2_)
749 goSparse = 1;
750 } else {
751 if (numberNonZero < sparseThreshold_)
752 goSparse = 2;
753 }
754 }
755 #ifndef COIN_FAST_CODE
756 assert(slackValue_ == -1.0);
757 #endif
758 if (!goSparse && numberRows_ < 1000) {
759 double *COIN_RESTRICT arrayFT = regionFT->denseVector();
760 int *COIN_RESTRICT indexFT = regionFT->getIndices();
761 int numberNonZeroFT;
762 double *COIN_RESTRICT arrayUpdate = regionUpdate->denseVector();
763 int *COIN_RESTRICT indexUpdate = regionUpdate->getIndices();
764 int numberNonZeroUpdate;
765 updateTwoColumnsUDensish(numberNonZeroFT, arrayFT, indexFT,
766 numberNonZeroUpdate, arrayUpdate, indexUpdate);
767 regionFT->setNumElements(numberNonZeroFT);
768 regionUpdate->setNumElements(numberNonZeroUpdate);
769 if (collectStatistics_) {
770 ftranCountAfterU_ += numberNonZeroFT + numberNonZeroUpdate;
771 #ifdef CLP_FACTORIZATION_INSTRUMENT
772 scaledLengthDense += numberDense_ * numberNonZeroFT;
773 scaledLengthDenseSquared += numberDense_ * numberDense_ * numberNonZeroFT;
774 scaledLengthL += lengthL_ * numberNonZeroFT;
775 scaledLengthR += lengthR_ * numberNonZeroFT;
776 scaledLengthU += lengthU_ * numberNonZeroFT;
777 scaledLengthDense += numberDense_ * numberNonZeroUpdate;
778 scaledLengthDenseSquared += numberDense_ * numberDense_ * numberNonZeroUpdate;
779 scaledLengthL += lengthL_ * numberNonZeroUpdate;
780 scaledLengthR += lengthR_ * numberNonZeroUpdate;
781 scaledLengthU += lengthU_ * numberNonZeroUpdate;
782 #endif
783 }
784 } else {
785 // sparse
786 updateColumnU(regionFT, regionIndex);
787 updateColumnU(regionUpdate, regionUpdate->getIndices());
788 }
789 permuteBack(regionFT, regionSparse2);
790 if (!noPermuteRegion3) {
791 permuteBack(regionUpdate, regionSparse3);
792 }
793 #ifdef COIN_DEBUG
794 int n2 = regionSparse2->getNumElements();
795 regionSparse1->checkClean();
796 int n2a = updateColumnFT(regionSparse1, &save2);
797 assert(n2 == n2a);
798 {
799 int j;
800 double *regionA = save2.denseVector();
801 int *indexA = save2.getIndices();
802 double *regionB = regionSparse2->denseVector();
803 int *indexB = regionSparse2->getIndices();
804 for (j = 0; j < n2; j++) {
805 int k = indexA[j];
806 assert(k == indexB[j]);
807 CoinFactorizationDouble value = regionA[j];
808 assert(value == regionB[j]);
809 }
810 }
811 updateColumn(&save3,
812 &save3,
813 noPermuteRegion3);
814 int n3 = regionSparse3->getNumElements();
815 assert(n3 == save3.getNumElements());
816 {
817 int j;
818 double *regionA = save3.denseVector();
819 int *indexA = save3.getIndices();
820 double *regionB = regionSparse3->denseVector();
821 int *indexB = regionSparse3->getIndices();
822 for (j = 0; j < n3; j++) {
823 int k = indexA[j];
824 assert(k == indexB[j]);
825 CoinFactorizationDouble value = regionA[k];
826 assert(value == regionB[k]);
827 }
828 }
829 //*regionSparse2=save2;
830 //*regionSparse3=save3;
831 printf("REGION2 %d els\n", regionSparse2->getNumElements());
832 regionSparse2->print();
833 printf("REGION3 %d els\n", regionSparse3->getNumElements());
834 regionSparse3->print();
835 #endif
836 return regionSparse2->getNumElements();
837 #else
838 int returnCode = updateColumnFT(regionSparse1,
839 regionSparse2);
840 assert(noPermuteRegion3);
841 updateColumn(regionSparse3,
842 regionSparse3,
843 noPermuteRegion3);
844 #ifdef CLP_FACTORIZATION_INSTRUMENT
845 numberUpdateTwoFT++;
846 timeInUpdateTwoFT += CoinCpuTime() - startTimeX;
847 averageLengthR += 2 * lengthR_;
848 averageLengthU += 2 * lengthU_;
849 averageLengthL += 2 * lengthL_;
850 #endif
851 //printf("REGION2 %d els\n",regionSparse2->getNumElements());
852 //regionSparse2->print();
853 //printf("REGION3 %d els\n",regionSparse3->getNumElements());
854 //regionSparse3->print();
855 return returnCode;
856 #endif
857 }
858 // Updates part of 2 columns (FTRANU) real work
updateTwoColumnsUDensish(int & numberNonZero1,double * COIN_RESTRICT region1,int * COIN_RESTRICT index1,int & numberNonZero2,double * COIN_RESTRICT region2,int * COIN_RESTRICT index2) const859 void CoinFactorization::updateTwoColumnsUDensish(
860 int &numberNonZero1,
861 double *COIN_RESTRICT region1,
862 int *COIN_RESTRICT index1,
863 int &numberNonZero2,
864 double *COIN_RESTRICT region2,
865 int *COIN_RESTRICT index2) const
866 {
867 double tolerance = zeroTolerance_;
868 const int *COIN_RESTRICT startColumn = startColumnU_.array();
869 const int *COIN_RESTRICT indexRow = indexRowU_.array();
870 const CoinFactorizationDouble *COIN_RESTRICT element = elementU_.array();
871 int numberNonZeroA = 0;
872 int numberNonZeroB = 0;
873 const int *numberInColumn = numberInColumn_.array();
874 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
875
876 for (int i = numberU_ - 1; i >= numberSlacks_; i--) {
877 CoinFactorizationDouble pivotValue2 = region2[i];
878 region2[i] = 0.0;
879 CoinFactorizationDouble pivotValue1 = region1[i];
880 region1[i] = 0.0;
881 if (fabs(pivotValue2) > tolerance) {
882 int start = startColumn[i];
883 const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
884 const int *COIN_RESTRICT thisIndex = indexRow + start;
885 if (fabs(pivotValue1) <= tolerance) {
886 // just region 2
887 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
888 int iRow = thisIndex[j];
889 CoinFactorizationDouble value = thisElement[j];
890 #ifdef NO_LOAD
891 region2[iRow] -= value * pivotValue2;
892 #else
893 CoinFactorizationDouble regionValue2 = region2[iRow];
894 region2[iRow] = regionValue2 - value * pivotValue2;
895 #endif
896 }
897 pivotValue2 *= pivotRegion[i];
898 region2[i] = pivotValue2;
899 index2[numberNonZeroB++] = i;
900 } else {
901 // both
902 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
903 int iRow = thisIndex[j];
904 CoinFactorizationDouble value = thisElement[j];
905 #ifdef NO_LOAD
906 region1[iRow] -= value * pivotValue1;
907 region2[iRow] -= value * pivotValue2;
908 #else
909 CoinFactorizationDouble regionValue1 = region1[iRow];
910 CoinFactorizationDouble regionValue2 = region2[iRow];
911 region1[iRow] = regionValue1 - value * pivotValue1;
912 region2[iRow] = regionValue2 - value * pivotValue2;
913 #endif
914 }
915 pivotValue1 *= pivotRegion[i];
916 pivotValue2 *= pivotRegion[i];
917 region1[i] = pivotValue1;
918 index1[numberNonZeroA++] = i;
919 region2[i] = pivotValue2;
920 index2[numberNonZeroB++] = i;
921 }
922 } else if (fabs(pivotValue1) > tolerance) {
923 int start = startColumn[i];
924 const CoinFactorizationDouble *COIN_RESTRICT thisElement = element + start;
925 const int *COIN_RESTRICT thisIndex = indexRow + start;
926 // just region 1
927 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
928 int iRow = thisIndex[j];
929 CoinFactorizationDouble value = thisElement[j];
930 #ifdef NO_LOAD
931 region1[iRow] -= value * pivotValue1;
932 #else
933 CoinFactorizationDouble regionValue1 = region1[iRow];
934 region1[iRow] = regionValue1 - value * pivotValue1;
935 #endif
936 }
937 pivotValue1 *= pivotRegion[i];
938 region1[i] = pivotValue1;
939 index1[numberNonZeroA++] = i;
940 }
941 }
942 // Slacks
943
944 for (int i = numberSlacks_ - 1; i >= 0; i--) {
945 double value2 = region2[i];
946 double value1 = region1[i];
947 bool value1NonZero = (value1 != 0.0);
948 if (fabs(value2) > tolerance) {
949 region2[i] = -value2;
950 index2[numberNonZeroB++] = i;
951 } else {
952 region2[i] = 0.0;
953 }
954 if (value1NonZero) {
955 index1[numberNonZeroA] = i;
956 if (fabs(value1) > tolerance) {
957 region1[i] = -value1;
958 numberNonZeroA++;
959 } else {
960 region1[i] = 0.0;
961 }
962 }
963 }
964 numberNonZero1 = numberNonZeroA;
965 numberNonZero2 = numberNonZeroB;
966 }
967 #ifdef COIN_FACTORIZATION_DIAGNOSE
968 static int numberTimesX = 0;
969 static int numberSparseX = 0;
970 double sumNumberSparseX = 0.0;
971 static int numberSparsishX = 0;
972 double sumNumberSparsishX = 0.0;
973 static int numberDensishX = 0;
974 double sumNumberDensishX = 0.0;
975 #endif
976 // updateColumnU. Updates part of column (FTRANU)
updateColumnU(CoinIndexedVector * regionSparse,int * indexIn) const977 void CoinFactorization::updateColumnU(CoinIndexedVector *regionSparse,
978 int *indexIn) const
979 {
980 int numberNonZero = regionSparse->getNumElements();
981
982 int goSparse;
983 // Guess at number at end
984 if (sparseThreshold_ > 0) {
985 if (ftranAverageAfterR_) {
986 int newNumber = static_cast< int >(numberNonZero * ftranAverageAfterU_);
987 if (newNumber < sparseThreshold_)
988 goSparse = 2;
989 else if (newNumber < sparseThreshold2_)
990 goSparse = 1;
991 else
992 goSparse = 0;
993 } else {
994 if (numberNonZero < sparseThreshold_)
995 goSparse = 2;
996 else
997 goSparse = 0;
998 }
999 } else {
1000 goSparse = 0;
1001 }
1002 #ifdef COIN_FACTORIZATION_DIAGNOSE
1003 numberTimesX++;
1004 if (!goSparse) {
1005 numberDensishX++;
1006 sumNumberDensishX += numberNonZero;
1007 } else if (goSparse == 1) {
1008 numberSparsishX++;
1009 sumNumberSparsishX += numberNonZero;
1010 } else {
1011 numberSparseX++;
1012 sumNumberSparseX += numberNonZero;
1013 }
1014 if ((numberTimesX % 1000) == 0) {
1015 double averageDensish = (numberDensishX) ? sumNumberDensishX / numberDensishX : 0.0;
1016 double averageSparsish = (numberSparsishX) ? sumNumberSparsishX / numberSparsishX : 0.0;
1017 double averageSparse = (numberSparseX) ? sumNumberSparseX / numberSparseX : 0.0;
1018 printf("sparsity D,ish,S (%d,%g) , (%d,%g) , (%d,%g) - ftranFactor %g\n",
1019 numberDensishX, averageDensish,
1020 numberSparsishX, averageSparsish,
1021 numberSparseX, averageSparse, ftranAverageAfterU_);
1022 }
1023 #endif
1024 switch (goSparse) {
1025 case 0: // densish
1026 {
1027 double *region = regionSparse->denseVector();
1028 int *regionIndex = regionSparse->getIndices();
1029 int numberNonZero = updateColumnUDensish(region, regionIndex);
1030 regionSparse->setNumElements(numberNonZero);
1031 } break;
1032 case 1: // middling
1033 updateColumnUSparsish(regionSparse, indexIn);
1034 break;
1035 case 2: // sparse
1036 updateColumnUSparse(regionSparse, indexIn);
1037 break;
1038 }
1039 if (collectStatistics_) {
1040 ftranCountAfterU_ += regionSparse->getNumElements();
1041 #ifdef CLP_FACTORIZATION_INSTRUMENT
1042 int numberNonZero = regionSparse->getNumElements();
1043 scaledLengthDense += numberDense_ * numberNonZero;
1044 scaledLengthDenseSquared += numberDense_ * numberDense_ * numberNonZero;
1045 scaledLengthL += lengthL_ * numberNonZero;
1046 scaledLengthR += lengthR_ * numberNonZero;
1047 scaledLengthU += lengthU_ * numberNonZero;
1048 #endif
1049 }
1050 }
1051 #ifdef COIN_DEVELOP
1052 double ncall_DZ = 0.0;
1053 double nrow_DZ = 0.0;
1054 double nslack_DZ = 0.0;
1055 double nU_DZ = 0.0;
1056 double nnz_DZ = 0.0;
1057 double nDone_DZ = 0.0;
1058 #endif
1059 // Updates part of column (FTRANU) real work
updateColumnUDensish(double * COIN_RESTRICT region,int * COIN_RESTRICT regionIndex) const1060 int CoinFactorization::updateColumnUDensish(double *COIN_RESTRICT region,
1061 int *COIN_RESTRICT regionIndex) const
1062 {
1063 double tolerance = zeroTolerance_;
1064 const int *startColumn = startColumnU_.array();
1065 const int *indexRow = indexRowU_.array();
1066 const CoinFactorizationDouble *element = elementU_.array();
1067 int numberNonZero = 0;
1068 const int *numberInColumn = numberInColumn_.array();
1069 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
1070 #ifdef COIN_DEVELOP
1071 ncall_DZ++;
1072 nrow_DZ += numberRows_;
1073 nslack_DZ += numberSlacks_;
1074 nU_DZ += numberU_;
1075 #endif
1076
1077 for (int i = numberU_ - 1; i >= numberSlacks_; i--) {
1078 CoinFactorizationDouble pivotValue = region[i];
1079 if (pivotValue) {
1080 #ifdef COIN_DEVELOP
1081 nnz_DZ++;
1082 #endif
1083 region[i] = 0.0;
1084 if (fabs(pivotValue) > tolerance) {
1085 int start = startColumn[i];
1086 const CoinFactorizationDouble *thisElement = element + start;
1087 const int *thisIndex = indexRow + start;
1088 #ifdef COIN_DEVELOP
1089 nDone_DZ += numberInColumn[i];
1090 #endif
1091 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
1092 int iRow = thisIndex[j];
1093 CoinFactorizationDouble regionValue = region[iRow];
1094 CoinFactorizationDouble value = thisElement[j];
1095 region[iRow] = regionValue - value * pivotValue;
1096 }
1097 pivotValue *= pivotRegion[i];
1098 region[i] = pivotValue;
1099 regionIndex[numberNonZero++] = i;
1100 }
1101 }
1102 }
1103
1104 // now do slacks
1105 #ifndef COIN_FAST_CODE
1106 if (slackValue_ == -1.0) {
1107 #endif
1108 #if 0
1109 // Could skew loop to pick up next one earlier
1110 // might improve pipelining
1111 for (int i = numberSlacks_-1; i>2;i-=2) {
1112 double value0 = region[i];
1113 double absValue0 = fabs ( value0 );
1114 double value1 = region[i-1];
1115 double absValue1 = fabs ( value1 );
1116 if ( value0 ) {
1117 if ( absValue0 > tolerance ) {
1118 region[i]=-value0;
1119 regionIndex[numberNonZero++]=i;
1120 } else {
1121 region[i]=0.0;
1122 }
1123 }
1124 if ( value1 ) {
1125 if ( absValue1 > tolerance ) {
1126 region[i-1]=-value1;
1127 regionIndex[numberNonZero++]=i-1;
1128 } else {
1129 region[i-1]=0.0;
1130 }
1131 }
1132 }
1133 for ( ; i>=0;i--) {
1134 double value = region[i];
1135 double absValue = fabs ( value );
1136 if ( value ) {
1137 if ( absValue > tolerance ) {
1138 region[i]=-value;
1139 regionIndex[numberNonZero++]=i;
1140 } else {
1141 region[i]=0.0;
1142 }
1143 }
1144 }
1145 #else
1146 for (int i = numberSlacks_ - 1; i >= 0; i--) {
1147 double value = region[i];
1148 if (value) {
1149 region[i] = -value;
1150 regionIndex[numberNonZero] = i;
1151 if (fabs(value) > tolerance)
1152 numberNonZero++;
1153 else
1154 region[i] = 0.0;
1155 }
1156 }
1157 #endif
1158 #ifndef COIN_FAST_CODE
1159 } else {
1160 assert(slackValue_ == 1.0);
1161 for (int i = numberSlacks_ - 1; i >= 0; i--) {
1162 double value = region[i];
1163 double absValue = fabs(value);
1164 if (value) {
1165 region[i] = 0.0;
1166 if (absValue > tolerance) {
1167 region[i] = value;
1168 regionIndex[numberNonZero++] = i;
1169 }
1170 }
1171 }
1172 }
1173 #endif
1174 return numberNonZero;
1175 }
1176 // updateColumnU. Updates part of column (FTRANU)
1177 /*
1178 Since everything is in order I should be able to do a better job of
1179 marking stuff - think. Also as L is static maybe I can do something
1180 better there (I know I could if I marked the depth of every element
1181 but that would lead to other inefficiencies.
1182 */
updateColumnUSparse(CoinIndexedVector * regionSparse,int * COIN_RESTRICT indexIn) const1183 void CoinFactorization::updateColumnUSparse(CoinIndexedVector *regionSparse,
1184 int *COIN_RESTRICT indexIn) const
1185 {
1186 int numberNonZero = regionSparse->getNumElements();
1187 int *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1188 double *COIN_RESTRICT region = regionSparse->denseVector();
1189 double tolerance = zeroTolerance_;
1190 const int *startColumn = startColumnU_.array();
1191 const int *indexRow = indexRowU_.array();
1192 const CoinFactorizationDouble *element = elementU_.array();
1193 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
1194 // use sparse_ as temporary area
1195 // mark known to be zero
1196 #if ABOCA_LITE_FACTORIZATION == 0
1197 #define sparseOffset 0
1198 #else
1199 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
1200 assert(!sparseOffset);
1201 #endif
1202 int *COIN_RESTRICT stack = sparse_.array() + sparseOffset; /* pivot */
1203 int *COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
1204 int *COIN_RESTRICT next = reinterpret_cast< int * >(list + maximumRowsExtra_); /* jnext */
1205 char *COIN_RESTRICT mark = reinterpret_cast< char * >(next + maximumRowsExtra_);
1206 #ifdef COIN_DEBUG
1207 for (int i = 0; i < maximumRowsExtra_; i++) {
1208 assert(!mark[i]);
1209 }
1210 #endif
1211
1212 // move slacks to end of stack list
1213 int *COIN_RESTRICT putLast = stack + maximumRowsExtra_;
1214 int *COIN_RESTRICT put = putLast;
1215
1216 const int *numberInColumn = numberInColumn_.array();
1217 int nList = 0;
1218 for (int i = 0; i < numberNonZero; i++) {
1219 int kPivot = indexIn[i];
1220 stack[0] = kPivot;
1221 int j = startColumn[kPivot] + numberInColumn[kPivot] - 1;
1222 int nStack = 1;
1223 next[0] = j;
1224 while (nStack) {
1225 /* take off stack */
1226 int kPivot = stack[--nStack];
1227 if (mark[kPivot] != 1) {
1228 j = next[nStack];
1229 if (j >= startColumn[kPivot]) {
1230 kPivot = indexRow[j--];
1231 /* put back on stack */
1232 next[nStack++] = j;
1233 if (!mark[kPivot]) {
1234 /* and new one */
1235 int numberIn = numberInColumn[kPivot];
1236 if (numberIn) {
1237 j = startColumn[kPivot] + numberIn - 1;
1238 stack[nStack] = kPivot;
1239 mark[kPivot] = 2;
1240 next[nStack++] = j;
1241 } else {
1242 // can do immediately
1243 /* finished so mark */
1244 mark[kPivot] = 1;
1245 if (kPivot >= numberSlacks_) {
1246 list[nList++] = kPivot;
1247 } else {
1248 // slack - put at end
1249 --put;
1250 *put = kPivot;
1251 }
1252 }
1253 }
1254 } else {
1255 /* finished so mark */
1256 mark[kPivot] = 1;
1257 if (kPivot >= numberSlacks_) {
1258 list[nList++] = kPivot;
1259 } else {
1260 // slack - put at end
1261 assert(!numberInColumn[kPivot]);
1262 --put;
1263 *put = kPivot;
1264 }
1265 }
1266 }
1267 }
1268 }
1269 #if 0
1270 {
1271 std::sort(list,list+nList);
1272 int i;
1273 int last;
1274 last =-1;
1275 for (i=0;i<nList;i++) {
1276 int k = list[i];
1277 assert (k>last);
1278 last=k;
1279 }
1280 std::sort(put,putLast);
1281 int n = putLast-put;
1282 last =-1;
1283 for (i=0;i<n;i++) {
1284 int k = put[i];
1285 assert (k>last);
1286 last=k;
1287 }
1288 }
1289 #endif
1290 numberNonZero = 0;
1291 for (int i = nList - 1; i >= 0; i--) {
1292 int iPivot = list[i];
1293 mark[iPivot] = 0;
1294 CoinFactorizationDouble pivotValue = region[iPivot];
1295 region[iPivot] = 0.0;
1296 if (fabs(pivotValue) > tolerance) {
1297 int start = startColumn[iPivot];
1298 int number = numberInColumn[iPivot];
1299
1300 int j;
1301 for (j = start; j < start + number; j++) {
1302 CoinFactorizationDouble value = element[j];
1303 int iRow = indexRow[j];
1304 region[iRow] -= value * pivotValue;
1305 }
1306 pivotValue *= pivotRegion[iPivot];
1307 region[iPivot] = pivotValue;
1308 regionIndex[numberNonZero++] = iPivot;
1309 }
1310 }
1311 // slacks
1312 #ifndef COIN_FAST_CODE
1313 if (slackValue_ == 1.0) {
1314 for (; put < putLast; put++) {
1315 int iPivot = *put;
1316 mark[iPivot] = 0;
1317 CoinFactorizationDouble pivotValue = region[iPivot];
1318 region[iPivot] = 0.0;
1319 if (fabs(pivotValue) > tolerance) {
1320 region[iPivot] = pivotValue;
1321 regionIndex[numberNonZero++] = iPivot;
1322 }
1323 }
1324 } else {
1325 #endif
1326 for (; put < putLast; put++) {
1327 int iPivot = *put;
1328 mark[iPivot] = 0;
1329 CoinFactorizationDouble pivotValue = region[iPivot];
1330 region[iPivot] = 0.0;
1331 if (fabs(pivotValue) > tolerance) {
1332 region[iPivot] = -pivotValue;
1333 regionIndex[numberNonZero++] = iPivot;
1334 }
1335 }
1336 #ifndef COIN_FAST_CODE
1337 }
1338 #endif
1339 regionSparse->setNumElements(numberNonZero);
1340 }
1341 // updateColumnU. Updates part of column (FTRANU)
1342 /*
1343 Since everything is in order I should be able to do a better job of
1344 marking stuff - think. Also as L is static maybe I can do something
1345 better there (I know I could if I marked the depth of every element
1346 but that would lead to other inefficiencies.
1347 */
1348 #ifdef COIN_DEVELOP
1349 double ncall_SZ = 0.0;
1350 double nrow_SZ = 0.0;
1351 double nslack_SZ = 0.0;
1352 double nU_SZ = 0.0;
1353 double nnz_SZ = 0.0;
1354 double nDone_SZ = 0.0;
1355 #endif
updateColumnUSparsish(CoinIndexedVector * regionSparse,int * COIN_RESTRICT indexIn) const1356 void CoinFactorization::updateColumnUSparsish(CoinIndexedVector *regionSparse,
1357 int *COIN_RESTRICT indexIn) const
1358 {
1359 int *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1360 // mark known to be zero
1361 #if ABOCA_LITE_FACTORIZATION == 0
1362 #define sparseOffset 0
1363 #else
1364 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
1365 assert(!sparseOffset);
1366 #endif
1367 int *COIN_RESTRICT stack = sparse_.array() + sparseOffset; /* pivot */
1368 int *COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
1369 int *COIN_RESTRICT next = reinterpret_cast< int * >(list + maximumRowsExtra_); /* jnext */
1370 CoinCheckZero *COIN_RESTRICT mark = reinterpret_cast< CoinCheckZero * >(next + maximumRowsExtra_);
1371 const int *numberInColumn = numberInColumn_.array();
1372 #ifdef COIN_DEBUG
1373 for (int i = 0; i < maximumRowsExtra_; i++) {
1374 assert(!mark[i]);
1375 }
1376 #endif
1377
1378 int nMarked = 0;
1379 int numberNonZero = regionSparse->getNumElements();
1380 double *COIN_RESTRICT region = regionSparse->denseVector();
1381 double tolerance = zeroTolerance_;
1382 const int *startColumn = startColumnU_.array();
1383 const int *indexRow = indexRowU_.array();
1384 const CoinFactorizationDouble *element = elementU_.array();
1385 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
1386 #ifdef COIN_DEVELOP
1387 ncall_SZ++;
1388 nrow_SZ += numberRows_;
1389 nslack_SZ += numberSlacks_;
1390 nU_SZ += numberU_;
1391 #endif
1392
1393 for (int ii = 0; ii < numberNonZero; ii++) {
1394 int iPivot = indexIn[ii];
1395 int iWord = iPivot >> CHECK_SHIFT;
1396 int iBit = iPivot - (iWord << CHECK_SHIFT);
1397 if (mark[iWord]) {
1398 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
1399 } else {
1400 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
1401 stack[nMarked++] = iWord;
1402 }
1403 }
1404 numberNonZero = 0;
1405 // First do down to convenient power of 2
1406 int jLast = (numberU_ - 1) >> CHECK_SHIFT;
1407 jLast = CoinMax((jLast << CHECK_SHIFT), static_cast< int >(numberSlacks_));
1408 int i;
1409 for (i = numberU_ - 1; i >= jLast; i--) {
1410 CoinFactorizationDouble pivotValue = region[i];
1411 region[i] = 0.0;
1412 if (fabs(pivotValue) > tolerance) {
1413 #ifdef COIN_DEVELOP
1414 nnz_SZ++;
1415 #endif
1416 int start = startColumn[i];
1417 const CoinFactorizationDouble *thisElement = element + start;
1418 const int *thisIndex = indexRow + start;
1419
1420 #ifdef COIN_DEVELOP
1421 nDone_SZ += numberInColumn[i];
1422 #endif
1423 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
1424 int iRow0 = thisIndex[j];
1425 CoinFactorizationDouble regionValue0 = region[iRow0];
1426 CoinFactorizationDouble value0 = thisElement[j];
1427 int iWord = iRow0 >> CHECK_SHIFT;
1428 int iBit = iRow0 - (iWord << CHECK_SHIFT);
1429 if (mark[iWord]) {
1430 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
1431 } else {
1432 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
1433 stack[nMarked++] = iWord;
1434 }
1435 region[iRow0] = regionValue0 - value0 * pivotValue;
1436 }
1437 pivotValue *= pivotRegion[i];
1438 region[i] = pivotValue;
1439 regionIndex[numberNonZero++] = i;
1440 }
1441 }
1442 int kLast = (numberSlacks_ + BITS_PER_CHECK - 1) >> CHECK_SHIFT;
1443 if (jLast > numberSlacks_) {
1444 // now do in chunks
1445 for (int k = (jLast >> CHECK_SHIFT) - 1; k >= kLast; k--) {
1446 unsigned int iMark = mark[k];
1447 if (iMark) {
1448 // something in chunk - do all (as imark may change)
1449 int iLast = k << CHECK_SHIFT;
1450 for (i = iLast + BITS_PER_CHECK - 1; i >= iLast; i--) {
1451 CoinFactorizationDouble pivotValue = region[i];
1452 if (pivotValue) {
1453 #ifdef COIN_DEVELOP
1454 nnz_SZ++;
1455 #endif
1456 region[i] = 0.0;
1457 if (fabs(pivotValue) > tolerance) {
1458 int start = startColumn[i];
1459 const CoinFactorizationDouble *thisElement = element + start;
1460 const int *thisIndex = indexRow + start;
1461 #ifdef COIN_DEVELOP
1462 nDone_SZ += numberInColumn[i];
1463 #endif
1464 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
1465 int iRow0 = thisIndex[j];
1466 CoinFactorizationDouble regionValue0 = region[iRow0];
1467 CoinFactorizationDouble value0 = thisElement[j];
1468 int iWord = iRow0 >> CHECK_SHIFT;
1469 int iBit = iRow0 - (iWord << CHECK_SHIFT);
1470 if (mark[iWord]) {
1471 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
1472 } else {
1473 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
1474 stack[nMarked++] = iWord;
1475 }
1476 region[iRow0] = regionValue0 - value0 * pivotValue;
1477 }
1478 pivotValue *= pivotRegion[i];
1479 region[i] = pivotValue;
1480 regionIndex[numberNonZero++] = i;
1481 }
1482 }
1483 }
1484 mark[k] = 0;
1485 }
1486 }
1487 i = (kLast << CHECK_SHIFT) - 1;
1488 }
1489 for (; i >= numberSlacks_; i--) {
1490 CoinFactorizationDouble pivotValue = region[i];
1491 region[i] = 0.0;
1492 if (fabs(pivotValue) > tolerance) {
1493 #ifdef COIN_DEVELOP
1494 nnz_SZ++;
1495 #endif
1496 int start = startColumn[i];
1497 const CoinFactorizationDouble *thisElement = element + start;
1498 const int *thisIndex = indexRow + start;
1499 #ifdef COIN_DEVELOP
1500 nDone_SZ += numberInColumn[i];
1501 #endif
1502 for (int j = numberInColumn[i] - 1; j >= 0; j--) {
1503 int iRow0 = thisIndex[j];
1504 CoinFactorizationDouble regionValue0 = region[iRow0];
1505 CoinFactorizationDouble value0 = thisElement[j];
1506 int iWord = iRow0 >> CHECK_SHIFT;
1507 int iBit = iRow0 - (iWord << CHECK_SHIFT);
1508 if (mark[iWord]) {
1509 mark[iWord] = static_cast< CoinCheckZero >(mark[iWord] | (1 << iBit));
1510 } else {
1511 mark[iWord] = static_cast< CoinCheckZero >(1 << iBit);
1512 stack[nMarked++] = iWord;
1513 }
1514 region[iRow0] = regionValue0 - value0 * pivotValue;
1515 }
1516 pivotValue *= pivotRegion[i];
1517 region[i] = pivotValue;
1518 regionIndex[numberNonZero++] = i;
1519 }
1520 }
1521
1522 if (numberSlacks_) {
1523 // now do slacks
1524 #ifndef COIN_FAST_CODE
1525 double factor = slackValue_;
1526 if (factor == 1.0) {
1527 // First do down to convenient power of 2
1528 int jLast = (numberSlacks_ - 1) >> CHECK_SHIFT;
1529 jLast = jLast << CHECK_SHIFT;
1530 for (i = numberSlacks_ - 1; i >= jLast; i--) {
1531 double value = region[i];
1532 double absValue = fabs(value);
1533 if (value) {
1534 region[i] = 0.0;
1535 if (absValue > tolerance) {
1536 region[i] = value;
1537 regionIndex[numberNonZero++] = i;
1538 }
1539 }
1540 }
1541 mark[jLast] = 0;
1542 // now do in chunks
1543 for (int k = (jLast >> CHECK_SHIFT) - 1; k >= 0; k--) {
1544 unsigned int iMark = mark[k];
1545 if (iMark) {
1546 // something in chunk - do all (as imark may change)
1547 int iLast = k << CHECK_SHIFT;
1548 i = iLast + BITS_PER_CHECK - 1;
1549 for (; i >= iLast; i--) {
1550 double value = region[i];
1551 double absValue = fabs(value);
1552 if (value) {
1553 region[i] = 0.0;
1554 if (absValue > tolerance) {
1555 region[i] = value;
1556 regionIndex[numberNonZero++] = i;
1557 }
1558 }
1559 }
1560 mark[k] = 0;
1561 }
1562 }
1563 } else {
1564 assert(factor == -1.0);
1565 #endif
1566 // First do down to convenient power of 2
1567 int jLast = (numberSlacks_ - 1) >> CHECK_SHIFT;
1568 jLast = jLast << CHECK_SHIFT;
1569 for (i = numberSlacks_ - 1; i >= jLast; i--) {
1570 double value = region[i];
1571 double absValue = fabs(value);
1572 if (value) {
1573 region[i] = 0.0;
1574 if (absValue > tolerance) {
1575 region[i] = -value;
1576 regionIndex[numberNonZero++] = i;
1577 }
1578 }
1579 }
1580 mark[jLast] = 0;
1581 // now do in chunks
1582 for (int k = (jLast >> CHECK_SHIFT) - 1; k >= 0; k--) {
1583 unsigned int iMark = mark[k];
1584 if (iMark) {
1585 // something in chunk - do all (as imark may change)
1586 int iLast = k << CHECK_SHIFT;
1587 i = iLast + BITS_PER_CHECK - 1;
1588 for (; i >= iLast; i--) {
1589 double value = region[i];
1590 double absValue = fabs(value);
1591 if (value) {
1592 region[i] = 0.0;
1593 if (absValue > tolerance) {
1594 region[i] = -value;
1595 regionIndex[numberNonZero++] = i;
1596 }
1597 }
1598 }
1599 mark[k] = 0;
1600 }
1601 }
1602 #ifndef COIN_FAST_CODE
1603 }
1604 #endif
1605 }
1606 regionSparse->setNumElements(numberNonZero);
1607 mark[(numberU_ - 1) >> CHECK_SHIFT] = 0;
1608 mark[numberSlacks_ >> CHECK_SHIFT] = 0;
1609 if (numberSlacks_)
1610 mark[(numberSlacks_ - 1) >> CHECK_SHIFT] = 0;
1611 #ifdef COIN_DEBUG
1612 for (i = 0; i < maximumRowsExtra_; i++) {
1613 assert(!mark[i]);
1614 }
1615 #endif
1616 }
1617 // updateColumnR. Updates part of column (FTRANR)
updateColumnR(CoinIndexedVector * regionSparse) const1618 void CoinFactorization::updateColumnR(CoinIndexedVector *regionSparse) const
1619 {
1620 double *COIN_RESTRICT region = regionSparse->denseVector();
1621 int *COIN_RESTRICT regionIndex = regionSparse->getIndices();
1622 int numberNonZero = regionSparse->getNumElements();
1623
1624 if (!numberR_)
1625 return; //return if nothing to do
1626 double tolerance = zeroTolerance_;
1627
1628 const int *startColumn = startColumnR_.array() - numberRows_;
1629 const int *indexRow = indexRowR_;
1630 const CoinFactorizationDouble *element = elementR_;
1631 const int *permute = permute_.array();
1632
1633 // Work out very dubious idea of what would be fastest
1634 int method = -1;
1635 // Size of R
1636 double sizeR = startColumnR_.array()[numberR_];
1637 // Average
1638 double averageR = sizeR / (static_cast< double >(numberRowsExtra_));
1639 // weights (relative to actual work)
1640 double setMark = 0.1; // setting mark
1641 double test1 = 1.0; // starting ftran (without testPivot)
1642 double testPivot = 2.0; // Seeing if zero etc
1643 double startDot = 2.0; // For starting dot product version
1644 // For final scan
1645 double final = numberNonZero * 1.0;
1646 double methodTime[3];
1647 // For second type
1648 methodTime[1] = numberPivots_ * (testPivot + ((static_cast< double >(numberNonZero)) / (static_cast< double >(numberRows_)) * averageR));
1649 methodTime[1] += numberNonZero * (test1 + averageR);
1650 // For first type
1651 methodTime[0] = methodTime[1] + (numberNonZero + numberPivots_) * setMark;
1652 methodTime[1] += numberNonZero * final;
1653 // third
1654 methodTime[2] = sizeR + numberPivots_ * startDot + numberNonZero * final;
1655 // switch off if necessary
1656 if (!numberInColumnPlus_.array()) {
1657 methodTime[0] = 1.0e100;
1658 methodTime[1] = 1.0e100;
1659 } else if (!sparse_.array()) {
1660 methodTime[0] = 1.0e100;
1661 }
1662 double best = 1.0e100;
1663 for (int i = 0; i < 3; i++) {
1664 if (methodTime[i] < best) {
1665 best = methodTime[i];
1666 method = i;
1667 }
1668 }
1669 assert(method >= 0);
1670 const int *numberInColumnPlus = numberInColumnPlus_.array();
1671 //if (method==1)
1672 //printf(" methods %g %g %g - chosen %d\n",methodTime[0],methodTime[1],methodTime[2],method);
1673
1674 switch (method) {
1675 case 0:
1676 #ifdef STACK
1677 {
1678 // use sparse_ as temporary area
1679 // mark known to be zero
1680 #if ABOCA_LITE_FACTORIZATION == 0
1681 #define sparseOffset 0
1682 #else
1683 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
1684 assert(!sparseOffset);
1685 #endif
1686 int *COIN_RESTRICT stack = sparse_.array() + sparseOffset; /* pivot */
1687 int *COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
1688 int *COIN_RESTRICT next = (int *)(list + maximumRowsExtra_); /* jnext */
1689 char *COIN_RESTRICT mark = (char *)(next + maximumRowsExtra_);
1690 // we have another copy of R in R
1691 const CoinFactorizationDouble *elementR = elementR_ + lengthAreaR_;
1692 const int *indexRowR = indexRowR_ + lengthAreaR_;
1693 const int *startR = startColumnR_.array() + maximumPivots_ + 1;
1694 int nList = 0;
1695 const int *permuteBack = permuteBack_.array();
1696 for (int k = 0; k < numberNonZero; k++) {
1697 int kPivot = regionIndex[k];
1698 if (!mark[kPivot]) {
1699 stack[0] = kPivot;
1700 int j = -10;
1701 next[0] = j;
1702 int nStack = 0;
1703 while (nStack >= 0) {
1704 /* take off stack */
1705 if (j >= startR[kPivot]) {
1706 int jPivot = indexRowR[j--];
1707 /* put back on stack */
1708 next[nStack] = j;
1709 if (!mark[jPivot]) {
1710 /* and new one */
1711 kPivot = jPivot;
1712 j = -10;
1713 stack[++nStack] = kPivot;
1714 mark[kPivot] = 1;
1715 next[nStack] = j;
1716 }
1717 } else if (j == -10) {
1718 // before first - see if followon
1719 int jPivot = permuteBack[kPivot];
1720 if (jPivot < numberRows_) {
1721 // no
1722 j = startR[kPivot] + numberInColumnPlus[kPivot] - 1;
1723 next[nStack] = j;
1724 } else {
1725 // add to list
1726 if (!mark[jPivot]) {
1727 /* and new one */
1728 kPivot = jPivot;
1729 j = -10;
1730 stack[++nStack] = kPivot;
1731 mark[kPivot] = 1;
1732 next[nStack] = j;
1733 } else {
1734 j = startR[kPivot] + numberInColumnPlus[kPivot] - 1;
1735 next[nStack] = j;
1736 }
1737 }
1738 } else {
1739 // finished
1740 list[nList++] = kPivot;
1741 mark[kPivot] = 1;
1742 --nStack;
1743 if (nStack >= 0) {
1744 kPivot = stack[nStack];
1745 j = next[nStack];
1746 }
1747 }
1748 }
1749 }
1750 }
1751 numberNonZero = 0;
1752 for (int i = nList - 1; i >= 0; i--) {
1753 int iPivot = list[i];
1754 mark[iPivot] = 0;
1755 CoinFactorizationDouble pivotValue;
1756 if (iPivot < numberRows_) {
1757 pivotValue = region[iPivot];
1758 } else {
1759 int before = permute[iPivot];
1760 pivotValue = region[iPivot] + region[before];
1761 region[before] = 0.0;
1762 }
1763 if (fabs(pivotValue) > tolerance) {
1764 region[iPivot] = pivotValue;
1765 int start = startR[iPivot];
1766 int number = numberInColumnPlus[iPivot];
1767 int end = start + number;
1768 int j;
1769 for (j = start; j < end; j++) {
1770 int iRow = indexRowR[j];
1771 CoinFactorizationDouble value = elementR[j];
1772 region[iRow] -= value * pivotValue;
1773 }
1774 regionIndex[numberNonZero++] = iPivot;
1775 } else {
1776 region[iPivot] = 0.0;
1777 }
1778 }
1779 }
1780 #else
1781 {
1782
1783 // use sparse_ as temporary area
1784 // mark known to be zero
1785 #if ABOCA_LITE_FACTORIZATION == 0
1786 #define sparseOffset 0
1787 #else
1788 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
1789 assert(!sparseOffset);
1790 #endif
1791 int *COIN_RESTRICT stack = sparse_.array() + sparseOffset; /* pivot */
1792 int *COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
1793 int *COIN_RESTRICT next = reinterpret_cast< int * >(list + maximumRowsExtra_); /* jnext */
1794 char *COIN_RESTRICT mark = reinterpret_cast< char * >(next + maximumRowsExtra_);
1795 // mark all rows which will be permuted
1796 for (int i = numberRows_; i < numberRowsExtra_; i++) {
1797 int iRow = permute[i];
1798 mark[iRow] = 1;
1799 }
1800 // we have another copy of R in R
1801 const CoinFactorizationDouble *elementR = elementR_ + lengthAreaR_;
1802 const int *indexRowR = indexRowR_ + lengthAreaR_;
1803 const int *startR = startColumnR_.array() + maximumPivots_ + 1;
1804 // For current list order does not matter as
1805 // only affects end
1806 int newNumber = 0;
1807 for (int i = 0; i < numberNonZero; i++) {
1808 int iRow = regionIndex[i];
1809 assert(region[iRow]);
1810 if (!mark[iRow])
1811 regionIndex[newNumber++] = iRow;
1812 int number = numberInColumnPlus[iRow];
1813 if (number) {
1814 CoinFactorizationDouble pivotValue = region[iRow];
1815 int start = startR[iRow];
1816 int end = start + number;
1817 for (int j = start; j < end; j++) {
1818 CoinFactorizationDouble value = elementR[j];
1819 int jRow = indexRowR[j];
1820 region[jRow] -= pivotValue * value;
1821 }
1822 }
1823 }
1824 numberNonZero = newNumber;
1825 for (int i = numberRows_; i < numberRowsExtra_; i++) {
1826 //move using permute_ (stored in inverse fashion)
1827 int iRow = permute[i];
1828 CoinFactorizationDouble pivotValue = region[iRow] + region[i];
1829 //zero out pre-permuted
1830 region[iRow] = 0.0;
1831 if (fabs(pivotValue) > tolerance) {
1832 region[i] = pivotValue;
1833 if (!mark[i])
1834 regionIndex[numberNonZero++] = i;
1835 int number = numberInColumnPlus[i];
1836 int start = startR[i];
1837 int end = start + number;
1838 for (int j = start; j < end; j++) {
1839 CoinFactorizationDouble value = elementR[j];
1840 int jRow = indexRowR[j];
1841 region[jRow] -= pivotValue * value;
1842 }
1843 } else {
1844 region[i] = 0.0;
1845 }
1846 mark[iRow] = 0;
1847 }
1848 }
1849 #endif
1850 break;
1851 case 1: {
1852 // no sparse region
1853 // we have another copy of R in R
1854 const CoinFactorizationDouble *elementR = elementR_ + lengthAreaR_;
1855 const int *indexRowR = indexRowR_ + lengthAreaR_;
1856 const int *startR = startColumnR_.array() + maximumPivots_ + 1;
1857 // For current list order does not matter as
1858 // only affects end
1859 for (int i = 0; i < numberNonZero; i++) {
1860 int iRow = regionIndex[i];
1861 assert(region[iRow]);
1862 int number = numberInColumnPlus[iRow];
1863 if (number) {
1864 CoinFactorizationDouble pivotValue = region[iRow];
1865 int start = startR[iRow];
1866 int end = start + number;
1867 for (int j = start; j < end; j++) {
1868 CoinFactorizationDouble value = elementR[j];
1869 int jRow = indexRowR[j];
1870 region[jRow] -= pivotValue * value;
1871 }
1872 }
1873 }
1874 for (int i = numberRows_; i < numberRowsExtra_; i++) {
1875 //move using permute_ (stored in inverse fashion)
1876 int iRow = permute[i];
1877 CoinFactorizationDouble pivotValue = region[iRow] + region[i];
1878 //zero out pre-permuted
1879 region[iRow] = 0.0;
1880 if (fabs(pivotValue) > tolerance) {
1881 region[i] = pivotValue;
1882 regionIndex[numberNonZero++] = i;
1883 int number = numberInColumnPlus[i];
1884 int start = startR[i];
1885 int end = start + number;
1886 for (int j = start; j < end; j++) {
1887 CoinFactorizationDouble value = elementR[j];
1888 int jRow = indexRowR[j];
1889 region[jRow] -= pivotValue * value;
1890 }
1891 } else {
1892 region[i] = 0.0;
1893 }
1894 }
1895 } break;
1896 case 2: {
1897 int start = startColumn[numberRows_];
1898 for (int i = numberRows_; i < numberRowsExtra_; i++) {
1899 //move using permute_ (stored in inverse fashion)
1900 int end = startColumn[i + 1];
1901 int iRow = permute[i];
1902 CoinFactorizationDouble pivotValue = region[iRow];
1903 //zero out pre-permuted
1904 region[iRow] = 0.0;
1905
1906 for (int j = start; j < end; j++) {
1907 CoinFactorizationDouble value = element[j];
1908 int jRow = indexRow[j];
1909 value *= region[jRow];
1910 pivotValue -= value;
1911 }
1912 start = end;
1913 if (fabs(pivotValue) > tolerance) {
1914 region[i] = pivotValue;
1915 regionIndex[numberNonZero++] = i;
1916 } else {
1917 region[i] = 0.0;
1918 }
1919 }
1920 } break;
1921 }
1922 if (method) {
1923 // pack down
1924 int n = numberNonZero;
1925 numberNonZero = 0;
1926 for (int i = 0; i < n; i++) {
1927 int indexValue = regionIndex[i];
1928 double value = region[indexValue];
1929 if (value)
1930 regionIndex[numberNonZero++] = indexValue;
1931 }
1932 }
1933 //set counts
1934 regionSparse->setNumElements(numberNonZero);
1935 }
1936 // updateColumnR. Updates part of column (FTRANR)
updateColumnRFT(CoinIndexedVector * regionSparse,int * COIN_RESTRICT regionIndex)1937 void CoinFactorization::updateColumnRFT(CoinIndexedVector *regionSparse,
1938 int *COIN_RESTRICT regionIndex)
1939 {
1940 double *COIN_RESTRICT region = regionSparse->denseVector();
1941 //int *regionIndex = regionSparse->getIndices ( );
1942 int *COIN_RESTRICT startColumnU = startColumnU_.array();
1943 int numberNonZero = regionSparse->getNumElements();
1944
1945 if (numberR_) {
1946 double tolerance = zeroTolerance_;
1947
1948 const int *startColumn = startColumnR_.array() - numberRows_;
1949 const int *indexRow = indexRowR_;
1950 const CoinFactorizationDouble *element = elementR_;
1951 const int *permute = permute_.array();
1952
1953 // Work out very dubious idea of what would be fastest
1954 int method = -1;
1955 // Size of R
1956 double sizeR = startColumnR_.array()[numberR_];
1957 // Average
1958 double averageR = sizeR / (static_cast< double >(numberRowsExtra_));
1959 // weights (relative to actual work)
1960 double setMark = 0.1; // setting mark
1961 double test1 = 1.0; // starting ftran (without testPivot)
1962 double testPivot = 2.0; // Seeing if zero etc
1963 double startDot = 2.0; // For starting dot product version
1964 // For final scan
1965 double final = numberNonZero * 1.0;
1966 double methodTime[3];
1967 // For second type
1968 methodTime[1] = numberPivots_ * (testPivot + ((static_cast< double >(numberNonZero)) / (static_cast< double >(numberRows_)) * averageR));
1969 methodTime[1] += numberNonZero * (test1 + averageR);
1970 // For first type
1971 methodTime[0] = methodTime[1] + (numberNonZero + numberPivots_) * setMark;
1972 methodTime[1] += numberNonZero * final;
1973 // third
1974 methodTime[2] = sizeR + numberPivots_ * startDot + numberNonZero * final;
1975 // switch off if necessary
1976 if (!numberInColumnPlus_.array()) {
1977 methodTime[0] = 1.0e100;
1978 methodTime[1] = 1.0e100;
1979 } else if (!sparse_.array()) {
1980 methodTime[0] = 1.0e100;
1981 }
1982 const int *numberInColumnPlus = numberInColumnPlus_.array();
1983 int *numberInColumn = numberInColumn_.array();
1984 // adjust for final scan
1985 methodTime[1] += final;
1986 double best = 1.0e100;
1987 for (int i = 0; i < 3; i++) {
1988 if (methodTime[i] < best) {
1989 best = methodTime[i];
1990 method = i;
1991 }
1992 }
1993 assert(method >= 0);
1994
1995 switch (method) {
1996 case 0: {
1997 // use sparse_ as temporary area
1998 // mark known to be zero
1999 #if ABOCA_LITE_FACTORIZATION == 0
2000 #define sparseOffset 0
2001 #else
2002 int sparseOffset = ((regionSparse->capacity() & 0x80000000) != 0) ? sparseOffset_ : 0;
2003 assert(!sparseOffset);
2004 #endif
2005 int *COIN_RESTRICT stack = sparse_.array() + sparseOffset; /* pivot */
2006 int *COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
2007 int *COIN_RESTRICT next = reinterpret_cast< int * >(list + maximumRowsExtra_); /* jnext */
2008 char *COIN_RESTRICT mark = reinterpret_cast< char * >(next + maximumRowsExtra_);
2009 // mark all rows which will be permuted
2010 for (int i = numberRows_; i < numberRowsExtra_; i++) {
2011 int iRow = permute[i];
2012 mark[iRow] = 1;
2013 }
2014 // we have another copy of R in R
2015 const CoinFactorizationDouble *elementR = elementR_ + lengthAreaR_;
2016 const int *indexRowR = indexRowR_ + lengthAreaR_;
2017 const int *startR = startColumnR_.array() + maximumPivots_ + 1;
2018 //save in U
2019 //in at end
2020 int iColumn = numberColumnsExtra_;
2021
2022 startColumnU[iColumn] = startColumnU[maximumColumnsExtra_];
2023 int start = startColumnU[iColumn];
2024
2025 //int * putIndex = indexRowU_ + start;
2026 CoinFactorizationDouble *COIN_RESTRICT putElement = elementU_.array() + start;
2027 // For current list order does not matter as
2028 // only affects end
2029 int newNumber = 0;
2030 for (int i = 0; i < numberNonZero; i++) {
2031 int iRow = regionIndex[i];
2032 CoinFactorizationDouble pivotValue = region[iRow];
2033 assert(region[iRow]);
2034 if (!mark[iRow]) {
2035 //putIndex[newNumber]=iRow;
2036 putElement[newNumber] = pivotValue;
2037 ;
2038 regionIndex[newNumber++] = iRow;
2039 }
2040 int number = numberInColumnPlus[iRow];
2041 if (number) {
2042 int start = startR[iRow];
2043 int end = start + number;
2044 for (int j = start; j < end; j++) {
2045 CoinFactorizationDouble value = elementR[j];
2046 int jRow = indexRowR[j];
2047 region[jRow] -= pivotValue * value;
2048 }
2049 }
2050 }
2051 numberNonZero = newNumber;
2052 for (int i = numberRows_; i < numberRowsExtra_; i++) {
2053 //move using permute_ (stored in inverse fashion)
2054 int iRow = permute[i];
2055 CoinFactorizationDouble pivotValue = region[iRow] + region[i];
2056 //zero out pre-permuted
2057 region[iRow] = 0.0;
2058 if (fabs(pivotValue) > tolerance) {
2059 region[i] = pivotValue;
2060 if (!mark[i]) {
2061 //putIndex[numberNonZero]=i;
2062 putElement[numberNonZero] = pivotValue;
2063 ;
2064 regionIndex[numberNonZero++] = i;
2065 }
2066 int number = numberInColumnPlus[i];
2067 int start = startR[i];
2068 int end = start + number;
2069 for (int j = start; j < end; j++) {
2070 CoinFactorizationDouble value = elementR[j];
2071 int jRow = indexRowR[j];
2072 region[jRow] -= pivotValue * value;
2073 }
2074 } else {
2075 region[i] = 0.0;
2076 }
2077 mark[iRow] = 0;
2078 }
2079 numberInColumn[iColumn] = numberNonZero;
2080 startColumnU[maximumColumnsExtra_] = start + numberNonZero;
2081 } break;
2082 case 1: {
2083 // no sparse region
2084 // we have another copy of R in R
2085 const CoinFactorizationDouble *elementR = elementR_ + lengthAreaR_;
2086 const int *indexRowR = indexRowR_ + lengthAreaR_;
2087 const int *startR = startColumnR_.array() + maximumPivots_ + 1;
2088 // For current list order does not matter as
2089 // only affects end
2090 for (int i = 0; i < numberNonZero; i++) {
2091 int iRow = regionIndex[i];
2092 assert(region[iRow]);
2093 int number = numberInColumnPlus[iRow];
2094 if (number) {
2095 CoinFactorizationDouble pivotValue = region[iRow];
2096 int start = startR[iRow];
2097 int end = start + number;
2098 for (int j = start; j < end; j++) {
2099 CoinFactorizationDouble value = elementR[j];
2100 int jRow = indexRowR[j];
2101 region[jRow] -= pivotValue * value;
2102 }
2103 }
2104 }
2105 for (int i = numberRows_; i < numberRowsExtra_; i++) {
2106 //move using permute_ (stored in inverse fashion)
2107 int iRow = permute[i];
2108 CoinFactorizationDouble pivotValue = region[iRow] + region[i];
2109 //zero out pre-permuted
2110 region[iRow] = 0.0;
2111 if (fabs(pivotValue) > tolerance) {
2112 region[i] = pivotValue;
2113 regionIndex[numberNonZero++] = i;
2114 int number = numberInColumnPlus[i];
2115 int start = startR[i];
2116 int end = start + number;
2117 for (int j = start; j < end; j++) {
2118 CoinFactorizationDouble value = elementR[j];
2119 int jRow = indexRowR[j];
2120 region[jRow] -= pivotValue * value;
2121 }
2122 } else {
2123 region[i] = 0.0;
2124 }
2125 }
2126 } break;
2127 case 2: {
2128 int start = startColumn[numberRows_];
2129 for (int i = numberRows_; i < numberRowsExtra_; i++) {
2130 //move using permute_ (stored in inverse fashion)
2131 int end = startColumn[i + 1];
2132 int iRow = permute[i];
2133 CoinFactorizationDouble pivotValue = region[iRow];
2134 //zero out pre-permuted
2135 region[iRow] = 0.0;
2136
2137 for (int j = start; j < end; j++) {
2138 CoinFactorizationDouble value = element[j];
2139 int jRow = indexRow[j];
2140 value *= region[jRow];
2141 pivotValue -= value;
2142 }
2143 start = end;
2144 if (fabs(pivotValue) > tolerance) {
2145 region[i] = pivotValue;
2146 regionIndex[numberNonZero++] = i;
2147 } else {
2148 region[i] = 0.0;
2149 }
2150 }
2151 } break;
2152 }
2153 if (method) {
2154 // pack down
2155 int n = numberNonZero;
2156 numberNonZero = 0;
2157 //save in U
2158 //in at end
2159 int iColumn = numberColumnsExtra_;
2160
2161 assert(startColumnU[iColumn] == startColumnU[maximumColumnsExtra_]);
2162 int start = startColumnU[iColumn];
2163
2164 int *COIN_RESTRICT putIndex = indexRowU_.array() + start;
2165 CoinFactorizationDouble *COIN_RESTRICT putElement = elementU_.array() + start;
2166 for (int i = 0; i < n; i++) {
2167 int indexValue = regionIndex[i];
2168 double value = region[indexValue];
2169 if (value) {
2170 putIndex[numberNonZero] = indexValue;
2171 putElement[numberNonZero] = value;
2172 regionIndex[numberNonZero++] = indexValue;
2173 }
2174 }
2175 numberInColumn[iColumn] = numberNonZero;
2176 startColumnU[maximumColumnsExtra_] = start + numberNonZero;
2177 }
2178 //set counts
2179 regionSparse->setNumElements(numberNonZero);
2180 } else {
2181 // No R but we still need to save column
2182 //save in U
2183 //in at end
2184 int *COIN_RESTRICT numberInColumn = numberInColumn_.array();
2185 numberNonZero = regionSparse->getNumElements();
2186 int iColumn = numberColumnsExtra_;
2187
2188 assert(startColumnU[iColumn] == startColumnU[maximumColumnsExtra_]);
2189 int start = startColumnU[iColumn];
2190 numberInColumn[iColumn] = numberNonZero;
2191 startColumnU[maximumColumnsExtra_] = start + numberNonZero;
2192
2193 int *COIN_RESTRICT putIndex = indexRowU_.array() + start;
2194 CoinFactorizationDouble *COIN_RESTRICT putElement = elementU_.array() + start;
2195 for (int i = 0; i < numberNonZero; i++) {
2196 int indexValue = regionIndex[i];
2197 double value = region[indexValue];
2198 putIndex[i] = indexValue;
2199 putElement[i] = value;
2200 }
2201 }
2202 }
2203 /* Updates one column (FTRAN) from region2 and permutes.
2204 region1 starts as zero
2205 Note - if regionSparse2 packed on input - will be packed on output
2206 - returns un-permuted result in region2 and region1 is zero */
updateColumnFT(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2)2207 int CoinFactorization::updateColumnFT(CoinIndexedVector *regionSparse,
2208 CoinIndexedVector *regionSparse2)
2209 {
2210 #ifdef CLP_FACTORIZATION_INSTRUMENT
2211 double startTimeX = CoinCpuTime();
2212 #endif
2213 //permute and move indices into index array
2214 int *COIN_RESTRICT regionIndex = regionSparse->getIndices();
2215 int numberNonZero = regionSparse2->getNumElements();
2216 const int *permute = permute_.array();
2217 int *COIN_RESTRICT index = regionSparse2->getIndices();
2218 double *COIN_RESTRICT region = regionSparse->denseVector();
2219 double *COIN_RESTRICT array = regionSparse2->denseVector();
2220 int *COIN_RESTRICT startColumnU = startColumnU_.array();
2221 bool doFT = doForrestTomlin_;
2222 // see if room
2223 if (doFT) {
2224 int iColumn = numberColumnsExtra_;
2225
2226 startColumnU[iColumn] = startColumnU[maximumColumnsExtra_];
2227 int start = startColumnU[iColumn];
2228 int space = lengthAreaU_ - (start + numberRowsExtra_);
2229 doFT = space >= 0;
2230 if (doFT) {
2231 regionIndex = indexRowU_.array() + start;
2232 } else {
2233 startColumnU[maximumColumnsExtra_] = lengthAreaU_ + 1;
2234 }
2235 }
2236
2237 #ifndef CLP_FACTORIZATION
2238 bool packed = regionSparse2->packedMode();
2239 if (packed) {
2240 #else
2241 assert(regionSparse2->packedMode());
2242 #endif
2243 for (int j = 0; j < numberNonZero; j++) {
2244 int iRow = index[j];
2245 double value = array[j];
2246 array[j] = 0.0;
2247 iRow = permute[iRow];
2248 region[iRow] = value;
2249 regionIndex[j] = iRow;
2250 }
2251 #ifndef CLP_FACTORIZATION
2252 } else {
2253 for (int j = 0; j < numberNonZero; j++) {
2254 int iRow = index[j];
2255 double value = array[iRow];
2256 array[iRow] = 0.0;
2257 iRow = permute[iRow];
2258 region[iRow] = value;
2259 regionIndex[j] = iRow;
2260 }
2261 }
2262 #endif
2263 regionSparse->setNumElements(numberNonZero);
2264 if (collectStatistics_) {
2265 numberFtranCounts_++;
2266 ftranCountInput_ += numberNonZero;
2267 }
2268
2269 // ******* L
2270 #if 0
2271 {
2272 double *region = regionSparse->denseVector ( );
2273 //int *regionIndex = regionSparse->getIndices ( );
2274 int numberNonZero = regionSparse->getNumElements ( );
2275 for (int i=0;i<numberNonZero;i++) {
2276 int iRow = regionIndex[i];
2277 assert (region[iRow]);
2278 }
2279 }
2280 #endif
2281 updateColumnL(regionSparse, regionIndex);
2282 #if 0
2283 {
2284 double *region = regionSparse->denseVector ( );
2285 //int *regionIndex = regionSparse->getIndices ( );
2286 int numberNonZero = regionSparse->getNumElements ( );
2287 for (int i=0;i<numberNonZero;i++) {
2288 int iRow = regionIndex[i];
2289 assert (region[iRow]);
2290 }
2291 }
2292 #endif
2293 if (collectStatistics_)
2294 ftranCountAfterL_ += regionSparse->getNumElements();
2295 //permute extra
2296 //row bits here
2297 if (doFT)
2298 updateColumnRFT(regionSparse, regionIndex);
2299 else
2300 updateColumnR(regionSparse);
2301 if (collectStatistics_)
2302 ftranCountAfterR_ += regionSparse->getNumElements();
2303 // ******* U
2304 updateColumnU(regionSparse, regionIndex);
2305 if (!doForrestTomlin_) {
2306 // Do PFI after everything else
2307 updateColumnPFI(regionSparse);
2308 }
2309 permuteBack(regionSparse, regionSparse2);
2310 #ifdef CLP_FACTORIZATION_INSTRUMENT
2311 numberUpdateFT++;
2312 timeInUpdateFT += CoinCpuTime() - startTimeX;
2313 averageLengthR += lengthR_;
2314 averageLengthU += lengthU_;
2315 averageLengthL += lengthL_;
2316 #endif
2317 // will be negative if no room
2318 if (doFT)
2319 return regionSparse2->getNumElements();
2320 else
2321 return -regionSparse2->getNumElements();
2322 }
2323
2324 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2325 */
2326