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