1/* $Id: CoinFactorization4.cpp 1535 2012-06-14 07:15:54Z forrest $ */
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
40// getRowSpaceIterate. Gets space for one Row with given length
41//may have to do compression (returns true)
42//also moves existing vector
43bool
44CoinFactorization::getRowSpaceIterate ( int iRow,
45 int extraNeeded )
46{
47 const int * numberInRow = numberInRow_.array();
48 int number = numberInRow[iRow];
49 CoinBigIndex * COIN_RESTRICT startRow = startRowU_.array();
50 int * COIN_RESTRICT indexColumn = indexColumnU_.array();
51 CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnU_.array();
52 CoinBigIndex space = lengthAreaU_ - startRow[maximumRowsExtra_];
53 int * COIN_RESTRICT nextRow = nextRow_.array();
54 int * COIN_RESTRICT lastRow = lastRow_.array();
55 if ( space < extraNeeded + number + 2 ) {
56 //compression
57 int iRow = nextRow[maximumRowsExtra_];
58 CoinBigIndex put = 0;
59 while ( iRow != maximumRowsExtra_ ) {
60 //move
61 CoinBigIndex get = startRow[iRow];
62 CoinBigIndex getEnd = startRow[iRow] + numberInRow[iRow];
63
64 startRow[iRow] = put;
65 CoinBigIndex i;
66 for ( i = get; i < getEnd; i++ ) {
67 indexColumn[put] = indexColumn[i];
68 convertRowToColumn[put] = convertRowToColumn[i];
69 put++;
70 }
71 iRow = nextRow[iRow];
72 } /* endwhile */
73 numberCompressions_++;
74 startRow[maximumRowsExtra_] = put;
75 space = lengthAreaU_ - put;
76 if ( space < extraNeeded + number + 2 ) {
77 //need more space
78 //if we can allocate bigger then do so and copy
79 //if not then return so code can start again
80 status_ = -99;
81 return false; }
82 }
83 CoinBigIndex put = startRow[maximumRowsExtra_];
84 int next = nextRow[iRow];
85 int last = lastRow[iRow];
86
87 //out
88 nextRow[last] = next;
89 lastRow[next] = last;
90 //in at end
91 last = lastRow[maximumRowsExtra_];
92 nextRow[last] = iRow;
93 lastRow[maximumRowsExtra_] = iRow;
94 lastRow[iRow] = last;
95 nextRow[iRow] = maximumRowsExtra_;
96 //move
97 CoinBigIndex get = startRow[iRow];
98
99 int * indexColumnU = indexColumnU_.array();
100 startRow[iRow] = put;
101 while ( number ) {
102 number--;
103 indexColumnU[put] = indexColumnU[get];
104 convertRowToColumn[put] = convertRowToColumn[get];
105 put++;
106 get++;
107 } /* endwhile */
108 //add four for luck
109 startRow[maximumRowsExtra_] = put + extraNeeded + 4;
110 return true;
111}
112
113// getColumnSpaceIterateR. Gets space for one extra R element in Column
114//may have to do compression (returns true)
115//also moves existing vector
116bool
117CoinFactorization::getColumnSpaceIterateR ( int iColumn, double value,
118 int iRow)
119{
120 CoinFactorizationDouble * COIN_RESTRICT elementR = elementR_ + lengthAreaR_;
121 int * COIN_RESTRICT indexRowR = indexRowR_ + lengthAreaR_;
122 CoinBigIndex * COIN_RESTRICT startR = startColumnR_.array()+maximumPivots_+1;
123 int * COIN_RESTRICT numberInColumnPlus = numberInColumnPlus_.array();
124 int number = numberInColumnPlus[iColumn];
125 //*** modify so sees if can go in
126 //see if it can go in at end
127 int * COIN_RESTRICT nextColumn = nextColumn_.array();
128 int * COIN_RESTRICT lastColumn = lastColumn_.array();
129 if (lengthAreaR_-startR[maximumColumnsExtra_]<number+1) {
130 //compression
131 int jColumn = nextColumn[maximumColumnsExtra_];
132 CoinBigIndex put = 0;
133 while ( jColumn != maximumColumnsExtra_ ) {
134 //move
135 CoinBigIndex get;
136 CoinBigIndex getEnd;
137
138 get = startR[jColumn];
139 getEnd = get + numberInColumnPlus[jColumn];
140 startR[jColumn] = put;
141 CoinBigIndex i;
142 for ( i = get; i < getEnd; i++ ) {
143 indexRowR[put] = indexRowR[i];
144 elementR[put] = elementR[i];
145 put++;
146 }
147 jColumn = nextColumn[jColumn];
148 }
149 numberCompressions_++;
150 startR[maximumColumnsExtra_]=put;
151 }
152 // Still may not be room (as iColumn was still in)
153 if (lengthAreaR_-startR[maximumColumnsExtra_]<number+1)
154 return false;
155
156 int next = nextColumn[iColumn];
157 int last = lastColumn[iColumn];
158
159 //out
160 nextColumn[last] = next;
161 lastColumn[next] = last;
162
163 CoinBigIndex put = startR[maximumColumnsExtra_];
164 //in at end
165 last = lastColumn[maximumColumnsExtra_];
166 nextColumn[last] = iColumn;
167 lastColumn[maximumColumnsExtra_] = iColumn;
168 lastColumn[iColumn] = last;
169 nextColumn[iColumn] = maximumColumnsExtra_;
170
171 //move
172 CoinBigIndex get = startR[iColumn];
173 startR[iColumn] = put;
174 int i = 0;
175 for (i=0 ; i < number; i ++ ) {
176 elementR[put]= elementR[get];
177 indexRowR[put++] = indexRowR[get++];
178 }
179 //insert
180 elementR[put]=value;
181 indexRowR[put++]=iRow;
182 numberInColumnPlus[iColumn]++;
183 //add 4 for luck
184 startR[maximumColumnsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaR_);
185 return true;
186}
187int CoinFactorization::checkPivot(double saveFromU,
188 double oldPivot) const
189{
190 int status;
191 if ( fabs ( saveFromU ) > 1.0e-8 ) {
192 double checkTolerance;
193
194 if ( numberRowsExtra_ < numberRows_ + 2 ) {
195 checkTolerance = 1.0e-5;
196 } else if ( numberRowsExtra_ < numberRows_ + 10 ) {
197 checkTolerance = 1.0e-6;
198 } else if ( numberRowsExtra_ < numberRows_ + 50 ) {
199 checkTolerance = 1.0e-8;
200 } else {
201 checkTolerance = 1.0e-10;
202 }
203 checkTolerance *= relaxCheck_;
204 if ( fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < checkTolerance ) {
205 status = 0;
206 } else {
207#if COIN_DEBUG
208 std::cout <<"inaccurate pivot "<< oldPivot << " "
209 << saveFromU << std::endl;
210#endif
211 if ( fabs ( fabs ( oldPivot ) - fabs ( saveFromU ) ) < 1.0e-12 ||
212 fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < 1.0e-8 ) {
213 status = 1;
214 } else {
215 status = 2;
216 }
217 }
218 } else {
219 //error
220 status = 2;
221#if COIN_DEBUG
222 std::cout <<"inaccurate pivot "<< saveFromU / oldPivot
223 << " " << saveFromU << std::endl;
224#endif
225 }
226 return status;
227}
228// replaceColumn. Replaces one Column to basis
229// returns 0=OK, 1=Probably OK, 2=singular, 3=no room
230//partial update already in U
231int
232CoinFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
233 int pivotRow,
234 double pivotCheck ,
235 bool checkBeforeModifying,
236 double )
237{
238 assert (numberU_<=numberRowsExtra_);
239 CoinBigIndex * COIN_RESTRICT startColumnU = startColumnU_.array();
240 CoinBigIndex * COIN_RESTRICT startColumn;
241 int * COIN_RESTRICT indexRow;
242 CoinFactorizationDouble * COIN_RESTRICT element;
243
244 //return at once if too many iterations
245 if ( numberColumnsExtra_ >= maximumColumnsExtra_ ) {
246 return 5;
247 }
248 if ( lengthAreaU_ < startColumnU[maximumColumnsExtra_] ) {
249 return 3;
250 }
251
252 int * COIN_RESTRICT numberInRow = numberInRow_.array();
253 int * COIN_RESTRICT numberInColumn = numberInColumn_.array();
254 int * COIN_RESTRICT numberInColumnPlus = numberInColumnPlus_.array();
255 int realPivotRow;
256 realPivotRow = pivotColumn_.array()[pivotRow];
257 //zeroed out region
258 double * COIN_RESTRICT region = regionSparse->denseVector ( );
259
260 element = elementU_.array();
261 //take out old pivot column
262
263 // If we have done no pivots then always check before modification
264 if (!numberPivots_)
265 checkBeforeModifying=true;
266
267 totalElements_ -= numberInColumn[realPivotRow];
268 CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegion_.array();
269 CoinFactorizationDouble oldPivot = pivotRegion[realPivotRow];
270 // for accuracy check
271 pivotCheck = pivotCheck / oldPivot;
272#if COIN_DEBUG>1
273 int checkNumber=1000000;
274 //if (numberL_) checkNumber=-1;
275 if (numberR_>=checkNumber) {
276 printf("pivot row %d, check %g - alpha region:\n",
277 realPivotRow,pivotCheck);
278 /*int i;
279 for (i=0;i<numberRows_;i++) {
280 if (pivotRegion[i])
281 printf("%d %g\n",i,pivotRegion[i]);
282 }*/
283 }
284#endif
285 pivotRegion[realPivotRow] = 0.0;
286
287 CoinBigIndex saveEnd = startColumnU[realPivotRow]
288 + numberInColumn[realPivotRow];
289 // not necessary at present - but take no chances for future
290 numberInColumn[realPivotRow] = 0;
291 //get entries in row (pivot not stored)
292 int numberNonZero = 0;
293 int * COIN_RESTRICT indexColumn = indexColumnU_.array();
294 CoinBigIndex * COIN_RESTRICT convertRowToColumn = convertRowToColumnU_.array();
295 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
296 CoinBigIndex * COIN_RESTRICT startRow = startRowU_.array();
297 CoinBigIndex start=0;
298 CoinBigIndex end=0;
299#if COIN_DEBUG>1
300 if (numberR_>=checkNumber)
301 printf("Before btranu\n");
302#endif
303
304#if COIN_ONE_ETA_COPY
305 if (convertRowToColumn) {
306#endif
307 start = startRow[realPivotRow];
308 end = start + numberInRow[realPivotRow];
309
310 int smallestIndex=numberRowsExtra_;
311 if (!checkBeforeModifying) {
312 for (CoinBigIndex i = start; i < end ; i ++ ) {
313 int iColumn = indexColumn[i];
314 assert (iColumn<numberRowsExtra_);
315 smallestIndex = CoinMin(smallestIndex,iColumn);
316 CoinBigIndex j = convertRowToColumn[i];
317
318 region[iColumn] = element[j];
319#if COIN_DEBUG>1
320 if (numberR_>=checkNumber)
321 printf("%d %g\n",iColumn,region[iColumn]);
322#endif
323 element[j] = 0.0;
324 regionIndex[numberNonZero++] = iColumn;
325 }
326 } else {
327 for (CoinBigIndex i = start; i < end ; i ++ ) {
328 int iColumn = indexColumn[i];
329 smallestIndex = CoinMin(smallestIndex,iColumn);
330 CoinBigIndex j = convertRowToColumn[i];
331
332 region[iColumn] = element[j];
333#if COIN_DEBUG>1
334 if (numberR_>=checkNumber)
335 printf("%d %g\n",iColumn,region[iColumn]);
336#endif
337 regionIndex[numberNonZero++] = iColumn;
338 }
339 }
340 //do BTRAN - finding first one to use
341 regionSparse->setNumElements ( numberNonZero );
342 updateColumnTransposeU ( regionSparse, smallestIndex );
343#if COIN_ONE_ETA_COPY
344 } else {
345 // use R to save where elements are
346 CoinBigIndex * saveWhere = NULL;
347 if (checkBeforeModifying) {
348 if ( lengthR_ + maximumRowsExtra_ +1>= lengthAreaR_ ) {
349 //not enough room
350 return 3;
351 }
352 saveWhere = indexRowR_+lengthR_;
353 }
354 replaceColumnU(regionSparse,saveWhere,
355 realPivotRow);
356 }
357#endif
358 numberNonZero = regionSparse->getNumElements ( );
359 CoinFactorizationDouble saveFromU = 0.0;
360
361 CoinBigIndex startU = startColumnU[numberColumnsExtra_];
362 int * COIN_RESTRICT indexU = &indexRowU_.array()[startU];
363 CoinFactorizationDouble * COIN_RESTRICT elementU = &elementU_.array()[startU];
364
365
366 // Do accuracy test here if caller is paranoid
367 if (checkBeforeModifying) {
368 double tolerance = zeroTolerance_;
369 int number = numberInColumn[numberColumnsExtra_];
370
371 for (CoinBigIndex i = 0; i < number; i++ ) {
372 int iRow = indexU[i];
373 //if (numberCompressions_==99&&lengthU_==278)
374 //printf("row %d saveFromU %g element %g region %g\n",
375 // iRow,saveFromU,elementU[i],region[iRow]);
376 if ( fabs ( elementU[i] ) > tolerance ) {
377 if ( iRow != realPivotRow ) {
378 saveFromU -= elementU[i] * region[iRow];
379 } else {
380 saveFromU += elementU[i];
381 }
382 }
383 }
384 //check accuracy
385 int status = checkPivot(saveFromU,pivotCheck);
386 if (status) {
387 // restore some things
388 pivotRegion[realPivotRow] = oldPivot;
389 number = saveEnd-startColumnU[realPivotRow];
390 totalElements_ += number;
391 numberInColumn[realPivotRow]=number;
392 regionSparse->clear();
393 return status;
394#if COIN_ONE_ETA_COPY
395 } else if (convertRowToColumn) {
396#else
397 } else {
398#endif
399 // do what we would have done by now
400 for (CoinBigIndex i = start; i < end ; i ++ ) {
401 CoinBigIndex j = convertRowToColumn[i];
402 element[j] = 0.0;
403 }
404#if COIN_ONE_ETA_COPY
405 } else {
406 // delete elements
407 // used R to save where elements are
408 CoinBigIndex * saveWhere = indexRowR_+lengthR_;
409 CoinFactorizationDouble *element = elementU_.array();
410 int n = saveWhere[0];
411 for (int i=0;i<n;i++) {
412 CoinBigIndex where = saveWhere[i+1];
413 element[where]=0.0;
414 }
415 //printf("deleting els\n");
416#endif
417 }
418 }
419 // Now zero out column of U
420 //take out old pivot column
421 for (CoinBigIndex i = startColumnU[realPivotRow]; i < saveEnd ; i ++ ) {
422 element[i] = 0.0;
423 }
424 //zero out pivot Row (before or after?)
425 //add to R
426 startColumn = startColumnR_.array();
427 indexRow = indexRowR_;
428 element = elementR_;
429 CoinBigIndex l = lengthR_;
430 int number = numberR_;
431
432 startColumn[number] = l; //for luck and first time
433 number++;
434 startColumn[number] = l + numberNonZero;
435 numberR_ = number;
436 lengthR_ = l + numberNonZero;
437 totalElements_ += numberNonZero;
438 if ( lengthR_ >= lengthAreaR_ ) {
439 //not enough room
440 regionSparse->clear();
441 return 3;
442 }
443#if COIN_DEBUG>1
444 if (numberR_>=checkNumber)
445 printf("After btranu\n");
446#endif
447 for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
448 int iRow = regionIndex[i];
449#if COIN_DEBUG>1
450 if (numberR_>=checkNumber)
451 printf("%d %g\n",iRow,region[iRow]);
452#endif
453
454 indexRow[l] = iRow;
455 element[l] = region[iRow];
456 l++;
457 }
458 int * nextRow;
459 int * lastRow;
460 int next;
461 int last;
462#if COIN_ONE_ETA_COPY
463 if (convertRowToColumn) {
464#endif
465 //take out row
466 nextRow = nextRow_.array();
467 lastRow = lastRow_.array();
468 next = nextRow[realPivotRow];
469 last = lastRow[realPivotRow];
470
471 nextRow[last] = next;
472 lastRow[next] = last;
473 numberInRow[realPivotRow]=0;
474#if COIN_DEBUG
475 nextRow[realPivotRow] = 777777;
476 lastRow[realPivotRow] = 777777;
477#endif
478#if COIN_ONE_ETA_COPY
479 }
480#endif
481 //do permute
482 permute_.array()[numberRowsExtra_] = realPivotRow;
483 // and other way
484 permuteBack_.array()[realPivotRow] = numberRowsExtra_;
485 permuteBack_.array()[numberRowsExtra_] = -1;
486 //and for safety
487 permute_.array()[numberRowsExtra_ + 1] = 0;
488
489 pivotColumn_.array()[pivotRow] = numberRowsExtra_;
490 pivotColumnBack()[numberRowsExtra_] = pivotRow;
491 startColumn = startColumnU;
492 indexRow = indexRowU_.array();
493 element = elementU_.array();
494
495 numberU_++;
496 number = numberInColumn[numberColumnsExtra_];
497
498 totalElements_ += number;
499 lengthU_ += number;
500 if ( lengthU_ >= lengthAreaU_ ) {
501 //not enough room
502 regionSparse->clear();
503 return 3;
504 }
505
506 saveFromU = 0.0;
507
508 //put in pivot
509 //add row counts
510
511#if COIN_DEBUG>1
512 if (numberR_>=checkNumber)
513 printf("On U\n");
514#endif
515#if COIN_ONE_ETA_COPY
516 if (convertRowToColumn) {
517#endif
518 for (CoinBigIndex i = 0; i < number; i++ ) {
519 int iRow = indexU[i];
520#if COIN_DEBUG>1
521 if (numberR_>=checkNumber)
522 printf("%d %g\n",iRow,elementU[i]);
523#endif
524
525 //assert ( fabs ( elementU[i] ) > zeroTolerance_ );
526 if ( iRow != realPivotRow ) {
527 int next = nextRow[iRow];
528 int iNumberInRow = numberInRow[iRow];
529 CoinBigIndex space;
530 CoinBigIndex put = startRow[iRow] + iNumberInRow;
531
532 space = startRow[next] - put;
533 if ( space <= 0 ) {
534 getRowSpaceIterate ( iRow, iNumberInRow + 4 );
535 put = startRow[iRow] + iNumberInRow;
536 }
537 indexColumn[put] = numberColumnsExtra_;
538 convertRowToColumn[put] = i + startU;
539 numberInRow[iRow] = iNumberInRow + 1;
540 saveFromU = saveFromU - elementU[i] * region[iRow];
541 } else {
542 //zero out and save
543 saveFromU += elementU[i];
544 elementU[i] = 0.0;
545 }
546 }
547 //in at end
548 last = lastRow[maximumRowsExtra_];
549 nextRow[last] = numberRowsExtra_;
550 lastRow[maximumRowsExtra_] = numberRowsExtra_;
551 lastRow[numberRowsExtra_] = last;
552 nextRow[numberRowsExtra_] = maximumRowsExtra_;
553 startRow[numberRowsExtra_] = startRow[maximumRowsExtra_];
554 numberInRow[numberRowsExtra_] = 0;
555#if COIN_ONE_ETA_COPY
556 } else {
557 //abort();
558 for (CoinBigIndex i = 0; i < number; i++ ) {
559 int iRow = indexU[i];
560#if COIN_DEBUG>1
561 if (numberR_>=checkNumber)
562 printf("%d %g\n",iRow,elementU[i]);
563#endif
564
565 if ( fabs ( elementU[i] ) > tolerance ) {
566 if ( iRow != realPivotRow ) {
567 saveFromU = saveFromU - elementU[i] * region[iRow];
568 } else {
569 //zero out and save
570 saveFromU += elementU[i];
571 elementU[i] = 0.0;
572 }
573 } else {
574 elementU[i] = 0.0;
575 }
576 }
577 }
578#endif
579 //column in at beginning (as empty)
580 int * COIN_RESTRICT nextColumn = nextColumn_.array();
581 int * COIN_RESTRICT lastColumn = lastColumn_.array();
582 next = nextColumn[maximumColumnsExtra_];
583 lastColumn[next] = numberColumnsExtra_;
584 nextColumn[maximumColumnsExtra_] = numberColumnsExtra_;
585 nextColumn[numberColumnsExtra_] = next;
586 lastColumn[numberColumnsExtra_] = maximumColumnsExtra_;
587 //check accuracy - but not if already checked (optimization problem)
588 int status = (checkBeforeModifying) ? 0 : checkPivot(saveFromU,pivotCheck);
589
590 if (status!=2) {
591
592 CoinFactorizationDouble pivotValue = 1.0 / saveFromU;
593
594 pivotRegion[numberRowsExtra_] = pivotValue;
595 //modify by pivot
596 for (CoinBigIndex i = 0; i < number; i++ ) {
597 elementU[i] *= pivotValue;
598 }
599 maximumU_ = CoinMax(maximumU_,startU+number);
600 numberRowsExtra_++;
601 numberColumnsExtra_++;
602 numberGoodU_++;
603 numberPivots_++;
604 }
605 if ( numberRowsExtra_ > numberRows_ + 50 ) {
606 CoinBigIndex extra = factorElements_ >> 1;
607
608 if ( numberRowsExtra_ > numberRows_ + 100 + numberRows_ / 500 ) {
609 if ( extra < 2 * numberRows_ ) {
610 extra = 2 * numberRows_;
611 }
612 } else {
613 if ( extra < 5 * numberRows_ ) {
614 extra = 5 * numberRows_;
615 }
616 }
617 CoinBigIndex added = totalElements_ - factorElements_;
618
619 if ( added > extra && added > ( factorElements_ ) << 1 && !status
620 && 3*totalElements_ > 2*(lengthAreaU_+lengthAreaL_)) {
621 status = 3;
622 if ( messageLevel_ & 4 ) {
623 std::cout << "Factorization has "<< totalElements_
624 << ", basis had " << factorElements_ <<std::endl;
625 }
626 }
627 }
628 if (numberInColumnPlus&&status<2) {
629 // we are going to put another copy of R in R
630 CoinFactorizationDouble * COIN_RESTRICT elementR = elementR_ + lengthAreaR_;
631 int * COIN_RESTRICT indexRowR = indexRowR_ + lengthAreaR_;
632 CoinBigIndex * COIN_RESTRICT startR = startColumnR_.array()+maximumPivots_+1;
633 int pivotRow = numberRowsExtra_-1;
634 for (CoinBigIndex i = 0; i < numberNonZero; i++ ) {
635 int iRow = regionIndex[i];
636 assert (pivotRow>iRow);
637 next = nextColumn[iRow];
638 CoinBigIndex space;
639 if (next!=maximumColumnsExtra_)
640 space = startR[next]-startR[iRow];
641 else
642 space = lengthAreaR_-startR[iRow];
643 int numberInR = numberInColumnPlus[iRow];
644 if (space>numberInR) {
645 // there is space
646 CoinBigIndex put=startR[iRow]+numberInR;
647 numberInColumnPlus[iRow]=numberInR+1;
648 indexRowR[put]=pivotRow;
649 elementR[put]=region[iRow];
650 //add 4 for luck
651 if (next==maximumColumnsExtra_)
652 startR[maximumColumnsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaR_);
653 } else {
654 // no space - do we shuffle?
655 if (!getColumnSpaceIterateR(iRow,region[iRow],pivotRow)) {
656 // printf("Need more space for R\n");
657 numberInColumnPlus_.conditionalDelete();
658 regionSparse->clear();
659 break;
660 }
661 }
662 region[iRow]=0.0;
663 }
664 regionSparse->setNumElements(0);
665 } else {
666 regionSparse->clear();
667 }
668 return status;
669}
670
671// updateColumnTranspose. Updates one column transpose (BTRAN)
672int
673CoinFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
674 CoinIndexedVector * regionSparse2 )
675 const
676{
677 //zero region
678 regionSparse->clear ( );
679 double * COIN_RESTRICT region = regionSparse->denseVector ( );
680 double * COIN_RESTRICT vector = regionSparse2->denseVector();
681 int * COIN_RESTRICT index = regionSparse2->getIndices();
682 int numberNonZero = regionSparse2->getNumElements();
683 const int * pivotColumn = pivotColumn_.array();
684
685 //move indices into index array
686 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
687 bool packed = regionSparse2->packedMode();
688 if (packed) {
689 for (int i = 0; i < numberNonZero; i ++ ) {
690 int iRow = index[i];
691 double value = vector[i];
692 iRow=pivotColumn[iRow];
693 vector[i]=0.0;
694 region[iRow] = value;
695 regionIndex[i] = iRow;
696 }
697 } else {
698 for (int i = 0; i < numberNonZero; i ++ ) {
699 int iRow = index[i];
700 double value = vector[iRow];
701 vector[iRow]=0.0;
702 iRow=pivotColumn[iRow];
703 region[iRow] = value;
704 regionIndex[i] = iRow;
705 }
706 }
707 regionSparse->setNumElements ( numberNonZero );
708 if (collectStatistics_) {
709 numberBtranCounts_++;
710 btranCountInput_ += static_cast<double> (numberNonZero);
711 }
712 if (!doForrestTomlin_) {
713 // Do PFI before everything else
714 updateColumnTransposePFI(regionSparse);
715 numberNonZero = regionSparse->getNumElements();
716 }
717 // ******* U
718 // Apply pivot region - could be combined for speed
719 CoinFactorizationDouble * COIN_RESTRICT pivotRegion = pivotRegion_.array();
720
721 int smallestIndex=numberRowsExtra_;
722 for (int j = 0; j < numberNonZero; j++ ) {
723 int iRow = regionIndex[j];
724 smallestIndex = CoinMin(smallestIndex,iRow);
725 region[iRow] *= pivotRegion[iRow];
726 }
727 updateColumnTransposeU ( regionSparse,smallestIndex );
728 if (collectStatistics_)
729 btranCountAfterU_ += static_cast<double> (regionSparse->getNumElements());
730 //permute extra
731 //row bits here
732 updateColumnTransposeR ( regionSparse );
733 // ******* L
734 updateColumnTransposeL ( regionSparse );
735 numberNonZero = regionSparse->getNumElements ( );
736 if (collectStatistics_)
737 btranCountAfterL_ += static_cast<double> (numberNonZero);
738 const int * permuteBack = pivotColumnBack();
739 int number=0;
740 if (packed) {
741 for (int i=0;i<numberNonZero;i++) {
742 int iRow=regionIndex[i];
743 double value = region[iRow];
744 region[iRow]=0.0;
745 //if (fabs(value)>zeroTolerance_) {
746 iRow=permuteBack[iRow];
747 vector[number]=value;
748 index[number++]=iRow;
749 //}
750 }
751 } else {
752 for (int i=0;i<numberNonZero;i++) {
753 int iRow=regionIndex[i];
754 double value = region[iRow];
755 region[iRow]=0.0;
756 //if (fabs(value)>zeroTolerance_) {
757 iRow=permuteBack[iRow];
758 vector[iRow]=value;
759 index[number++]=iRow;
760 //}
761 }
762 }
763 regionSparse->setNumElements(0);
764 regionSparse2->setNumElements(number);
765#ifdef COIN_DEBUG
766 for (i=0;i<numberRowsExtra_;i++) {
767 assert (!region[i]);
768 }
769#endif
770 return number;
771}
772
773/* Updates part of column transpose (BTRANU) when densish,
774 assumes index is sorted i.e. region is correct */
775void
776CoinFactorization::updateColumnTransposeUDensish
777 ( CoinIndexedVector * regionSparse,
778 int smallestIndex) const
779{
780 double * COIN_RESTRICT region = regionSparse->denseVector ( );
781 int numberNonZero = regionSparse->getNumElements ( );
782 double tolerance = zeroTolerance_;
783
784 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
785
786 const CoinBigIndex *startRow = startRowU_.array();
787
788 const CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
789 const int *indexColumn = indexColumnU_.array();
790
791 const CoinFactorizationDouble * element = elementU_.array();
792 int last = numberU_;
793
794 const int *numberInRow = numberInRow_.array();
795 numberNonZero = 0;
796 for (int i=smallestIndex ; i < last; i++ ) {
797 CoinFactorizationDouble pivotValue = region[i];
798 if ( fabs ( pivotValue ) > tolerance ) {
799 CoinBigIndex start = startRow[i];
800 int numberIn = numberInRow[i];
801 CoinBigIndex end = start + numberIn;
802 for (CoinBigIndex j = start ; j < end; j ++ ) {
803 int iRow = indexColumn[j];
804 CoinBigIndex getElement = convertRowToColumn[j];
805 CoinFactorizationDouble value = element[getElement];
806 region[iRow] -= value * pivotValue;
807 }
808 regionIndex[numberNonZero++] = i;
809 } else {
810 region[i] = 0.0;
811 }
812 }
813 //set counts
814 regionSparse->setNumElements ( numberNonZero );
815}
816/* Updates part of column transpose (BTRANU) when sparsish,
817 assumes index is sorted i.e. region is correct */
818void
819CoinFactorization::updateColumnTransposeUSparsish
820 ( CoinIndexedVector * regionSparse,
821 int smallestIndex) const
822{
823 double * COIN_RESTRICT region = regionSparse->denseVector ( );
824 int numberNonZero = regionSparse->getNumElements ( );
825 double tolerance = zeroTolerance_;
826
827 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
828
829 const CoinBigIndex *startRow = startRowU_.array();
830
831 const CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
832 const int *indexColumn = indexColumnU_.array();
833
834 const CoinFactorizationDouble * element = elementU_.array();
835 int last = numberU_;
836
837 const int *numberInRow = numberInRow_.array();
838
839 // mark known to be zero
840 int nInBig = sizeof(CoinBigIndex)/sizeof(int);
841 CoinCheckZero * COIN_RESTRICT mark = reinterpret_cast<CoinCheckZero *> (sparse_.array()+(2+nInBig)*maximumRowsExtra_);
842
843 for (int i=0;i<numberNonZero;i++) {
844 int iPivot=regionIndex[i];
845 int iWord = iPivot>>CHECK_SHIFT;
846 int iBit = iPivot-(iWord<<CHECK_SHIFT);
847 if (mark[iWord]) {
848 mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
849 } else {
850 mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
851 }
852 }
853
854 numberNonZero = 0;
855 // Find convenient power of 2
856 smallestIndex = smallestIndex >> CHECK_SHIFT;
857 int kLast = last>>CHECK_SHIFT;
858 // do in chunks
859
860 for (int k=smallestIndex;k<kLast;k++) {
861 unsigned int iMark = mark[k];
862 if (iMark) {
863 // something in chunk - do all (as imark may change)
864 int i = k<<CHECK_SHIFT;
865 int iLast = i+BITS_PER_CHECK;
866 for ( ; i < iLast; i++ ) {
867 CoinFactorizationDouble pivotValue = region[i];
868 if ( fabs ( pivotValue ) > tolerance ) {
869 CoinBigIndex start = startRow[i];
870 int numberIn = numberInRow[i];
871 CoinBigIndex end = start + numberIn;
872 for (CoinBigIndex j = start ; j < end; j ++ ) {
873 int iRow = indexColumn[j];
874 CoinBigIndex getElement = convertRowToColumn[j];
875 CoinFactorizationDouble value = element[getElement];
876 int iWord = iRow>>CHECK_SHIFT;
877 int iBit = iRow-(iWord<<CHECK_SHIFT);
878 if (mark[iWord]) {
879 mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
880 } else {
881 mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
882 }
883 region[iRow] -= value * pivotValue;
884 }
885 regionIndex[numberNonZero++] = i;
886 } else {
887 region[i] = 0.0;
888 }
889 }
890 mark[k]=0;
891 }
892 }
893 mark[kLast]=0;
894 for (int i = kLast<<CHECK_SHIFT; i < last; i++ ) {
895 CoinFactorizationDouble pivotValue = region[i];
896 if ( fabs ( pivotValue ) > tolerance ) {
897 CoinBigIndex start = startRow[i];
898 int numberIn = numberInRow[i];
899 CoinBigIndex end = start + numberIn;
900 for (CoinBigIndex j = start ; j < end; j ++ ) {
901 int iRow = indexColumn[j];
902 CoinBigIndex getElement = convertRowToColumn[j];
903 CoinFactorizationDouble value = element[getElement];
904
905 region[iRow] -= value * pivotValue;
906 }
907 regionIndex[numberNonZero++] = i;
908 } else {
909 region[i] = 0.0;
910 }
911 }
912#ifdef COIN_DEBUG
913 for (int i=0;i<maximumRowsExtra_;i++) {
914 assert (!mark[i]);
915 }
916#endif
917 //set counts
918 regionSparse->setNumElements ( numberNonZero );
919}
920/* Updates part of column transpose (BTRANU) when sparse,
921 assumes index is sorted i.e. region is correct */
922void
923CoinFactorization::updateColumnTransposeUSparse (
924 CoinIndexedVector * regionSparse) const
925{
926 double * COIN_RESTRICT region = regionSparse->denseVector ( );
927 int numberNonZero = regionSparse->getNumElements ( );
928 double tolerance = zeroTolerance_;
929
930 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
931 const CoinBigIndex *startRow = startRowU_.array();
932
933 const CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
934 const int *indexColumn = indexColumnU_.array();
935
936 const CoinFactorizationDouble * element = elementU_.array();
937
938 const int *numberInRow = numberInRow_.array();
939
940 // use sparse_ as temporary area
941 // mark known to be zero
942 int * COIN_RESTRICT stack = sparse_.array(); /* pivot */
943 int * COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
944 CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_); /* jnext */
945 char * COIN_RESTRICT mark = reinterpret_cast<char *> (next + maximumRowsExtra_);
946 int nList;
947#ifdef COIN_DEBUG
948 for (int i=0;i<maximumRowsExtra_;i++) {
949 assert (!mark[i]);
950 }
951#endif
952#if 0
953 {
954 int i;
955 for (i=0;i<numberRowsExtra_;i++) {
956 CoinBigIndex krs = startRow[i];
957 CoinBigIndex kre = krs + numberInRow[i];
958 CoinBigIndex k;
959 for (k=krs;k<kre;k++)
960 assert (indexColumn[k]>i);
961 }
962 }
963#endif
964 nList=0;
965 for (int k=0;k<numberNonZero;k++) {
966 int kPivot=regionIndex[k];
967 stack[0]=kPivot;
968 CoinBigIndex j=startRow[kPivot]+numberInRow[kPivot]-1;
969 next[0]=j;
970 int nStack=1;
971 while (nStack) {
972 /* take off stack */
973 kPivot=stack[--nStack];
974 if (mark[kPivot]!=1) {
975 j=next[nStack];
976 if (j>=startRow[kPivot]) {
977 kPivot=indexColumn[j--];
978 /* put back on stack */
979 next[nStack++] =j;
980 if (!mark[kPivot]) {
981 /* and new one */
982 j=startRow[kPivot]+numberInRow[kPivot]-1;
983 stack[nStack]=kPivot;
984 mark[kPivot]=2;
985 next[nStack++]=j;
986 }
987 } else {
988 // finished
989 list[nList++]=kPivot;
990 mark[kPivot]=1;
991 }
992 }
993 }
994 }
995 numberNonZero=0;
996 for (int i=nList-1;i>=0;i--) {
997 int iPivot = list[i];
998 mark[iPivot]=0;
999 CoinFactorizationDouble pivotValue = region[iPivot];
1000 if ( fabs ( pivotValue ) > tolerance ) {
1001 CoinBigIndex start = startRow[iPivot];
1002 int numberIn = numberInRow[iPivot];
1003 CoinBigIndex end = start + numberIn;
1004 for (CoinBigIndex j=start ; j < end; j ++ ) {
1005 int iRow = indexColumn[j];
1006 CoinBigIndex getElement = convertRowToColumn[j];
1007 CoinFactorizationDouble value = element[getElement];
1008 region[iRow] -= value * pivotValue;
1009 }
1010 regionIndex[numberNonZero++] = iPivot;
1011 } else {
1012 region[iPivot] = 0.0;
1013 }
1014 }
1015 //set counts
1016 regionSparse->setNumElements ( numberNonZero );
1017}
1018// updateColumnTransposeU. Updates part of column transpose (BTRANU)
1019//assumes index is sorted i.e. region is correct
1020//does not sort by sign
1021void
1022CoinFactorization::updateColumnTransposeU ( CoinIndexedVector * regionSparse,
1023 int smallestIndex) const
1024{
1025#if COIN_ONE_ETA_COPY
1026 CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
1027 if (!convertRowToColumn) {
1028 //abort();
1029 updateColumnTransposeUByColumn(regionSparse,smallestIndex);
1030 return;
1031 }
1032#endif
1033 int number = regionSparse->getNumElements ( );
1034 int goSparse;
1035 // Guess at number at end
1036 if (sparseThreshold_>0) {
1037 if (btranAverageAfterU_) {
1038 int newNumber = static_cast<int> (number*btranAverageAfterU_);
1039 if (newNumber< sparseThreshold_)
1040 goSparse = 2;
1041 else if (newNumber< sparseThreshold2_)
1042 goSparse = 1;
1043 else
1044 goSparse = 0;
1045 } else {
1046 if (number<sparseThreshold_)
1047 goSparse = 2;
1048 else
1049 goSparse = 0;
1050 }
1051 } else {
1052 goSparse=0;
1053 }
1054 switch (goSparse) {
1055 case 0: // densish
1056 updateColumnTransposeUDensish(regionSparse,smallestIndex);
1057 break;
1058 case 1: // middling
1059 updateColumnTransposeUSparsish(regionSparse,smallestIndex);
1060 break;
1061 case 2: // sparse
1062 updateColumnTransposeUSparse(regionSparse);
1063 break;
1064 }
1065}
1066
1067/* updateColumnTransposeLDensish.
1068 Updates part of column transpose (BTRANL) dense by column */
1069void
1070CoinFactorization::updateColumnTransposeLDensish
1071 ( CoinIndexedVector * regionSparse ) const
1072{
1073 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1074 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
1075 int numberNonZero;
1076 double tolerance = zeroTolerance_;
1077 int base;
1078 int first = -1;
1079
1080 numberNonZero=0;
1081 //scan
1082 for (first=numberRows_-1;first>=0;first--) {
1083 if (region[first])
1084 break;
1085 }
1086 if ( first >= 0 ) {
1087 base = baseL_;
1088 const CoinBigIndex * COIN_RESTRICT startColumn = startColumnL_.array();
1089 const int * COIN_RESTRICT indexRow = indexRowL_.array();
1090 const CoinFactorizationDouble * COIN_RESTRICT element = elementL_.array();
1091 int last = baseL_ + numberL_;
1092
1093 if ( first >= last ) {
1094 first = last - 1;
1095 }
1096 for (int i = first ; i >= base; i-- ) {
1097 CoinBigIndex j;
1098 CoinFactorizationDouble pivotValue = region[i];
1099 for ( j= startColumn[i] ; j < startColumn[i+1]; j++ ) {
1100 int iRow = indexRow[j];
1101 CoinFactorizationDouble value = element[j];
1102 pivotValue -= value * region[iRow];
1103 }
1104 if ( fabs ( pivotValue ) > tolerance ) {
1105 region[i] = pivotValue;
1106 regionIndex[numberNonZero++] = i;
1107 } else {
1108 region[i] = 0.0;
1109 }
1110 }
1111 //may have stopped early
1112 if ( first < base ) {
1113 base = first + 1;
1114 }
1115 if (base > 5) {
1116 int i=base-1;
1117 CoinFactorizationDouble pivotValue=region[i];
1118 bool store = fabs(pivotValue) > tolerance;
1119 for (; i > 0; i-- ) {
1120 bool oldStore = store;
1121 CoinFactorizationDouble oldValue = pivotValue;
1122 pivotValue = region[i-1];
1123 store = fabs ( pivotValue ) > tolerance ;
1124 if (!oldStore) {
1125 region[i] = 0.0;
1126 } else {
1127 region[i] = oldValue;
1128 regionIndex[numberNonZero++] = i;
1129 }
1130 }
1131 if (store) {
1132 region[0] = pivotValue;
1133 regionIndex[numberNonZero++] = 0;
1134 } else {
1135 region[0] = 0.0;
1136 }
1137 } else {
1138 for (int i = base -1 ; i >= 0; i-- ) {
1139 CoinFactorizationDouble pivotValue = region[i];
1140 if ( fabs ( pivotValue ) > tolerance ) {
1141 region[i] = pivotValue;
1142 regionIndex[numberNonZero++] = i;
1143 } else {
1144 region[i] = 0.0;
1145 }
1146 }
1147 }
1148 }
1149 //set counts
1150 regionSparse->setNumElements ( numberNonZero );
1151}
1152/* updateColumnTransposeLByRow.
1153 Updates part of column transpose (BTRANL) densish but by row */
1154void
1155CoinFactorization::updateColumnTransposeLByRow
1156 ( CoinIndexedVector * regionSparse ) const
1157{
1158 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1159 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
1160 int numberNonZero;
1161 double tolerance = zeroTolerance_;
1162 int first = -1;
1163
1164 // use row copy of L
1165 const CoinFactorizationDouble * element = elementByRowL_.array();
1166 const CoinBigIndex * startRow = startRowL_.array();
1167 const int * column = indexColumnL_.array();
1168 for (first=numberRows_-1;first>=0;first--) {
1169 if (region[first])
1170 break;
1171 }
1172 numberNonZero=0;
1173 for (int i=first;i>=0;i--) {
1174 CoinFactorizationDouble pivotValue = region[i];
1175 if ( fabs ( pivotValue ) > tolerance ) {
1176 regionIndex[numberNonZero++] = i;
1177 CoinBigIndex j;
1178 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
1179 int iRow = column[j];
1180 CoinFactorizationDouble value = element[j];
1181 region[iRow] -= pivotValue*value;
1182 }
1183 } else {
1184 region[i] = 0.0;
1185 }
1186 }
1187 //set counts
1188 regionSparse->setNumElements ( numberNonZero );
1189}
1190// Updates part of column transpose (BTRANL) when sparsish by row
1191void
1192CoinFactorization::updateColumnTransposeLSparsish
1193 ( CoinIndexedVector * regionSparse ) const
1194{
1195 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1196 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
1197 int numberNonZero = regionSparse->getNumElements();
1198 double tolerance = zeroTolerance_;
1199
1200 // use row copy of L
1201 const CoinFactorizationDouble * element = elementByRowL_.array();
1202 const CoinBigIndex * startRow = startRowL_.array();
1203 const int * column = indexColumnL_.array();
1204 // mark known to be zero
1205 int nInBig = sizeof(CoinBigIndex)/sizeof(int);
1206 CoinCheckZero * COIN_RESTRICT mark = reinterpret_cast<CoinCheckZero *> (sparse_.array()+(2+nInBig)*maximumRowsExtra_);
1207 for (int i=0;i<numberNonZero;i++) {
1208 int iPivot=regionIndex[i];
1209 int iWord = iPivot>>CHECK_SHIFT;
1210 int iBit = iPivot-(iWord<<CHECK_SHIFT);
1211 if (mark[iWord]) {
1212 mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
1213 } else {
1214 mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
1215 }
1216 }
1217 numberNonZero = 0;
1218 // First do down to convenient power of 2
1219 int jLast = (numberRows_-1)>>CHECK_SHIFT;
1220 jLast = (jLast<<CHECK_SHIFT);
1221 for (int i=numberRows_-1;i>=jLast;i--) {
1222 CoinFactorizationDouble pivotValue = region[i];
1223 if ( fabs ( pivotValue ) > tolerance ) {
1224 regionIndex[numberNonZero++] = i;
1225 CoinBigIndex j;
1226 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
1227 int iRow = column[j];
1228 CoinFactorizationDouble value = element[j];
1229 int iWord = iRow>>CHECK_SHIFT;
1230 int iBit = iRow-(iWord<<CHECK_SHIFT);
1231 if (mark[iWord]) {
1232 mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
1233 } else {
1234 mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
1235 }
1236 region[iRow] -= pivotValue*value;
1237 }
1238 } else {
1239 region[i] = 0.0;
1240 }
1241 }
1242 // and in chunks
1243 jLast = jLast>>CHECK_SHIFT;
1244 mark[jLast]=0;
1245 for (int k=jLast-1;k>=0;k--) {
1246 unsigned int iMark = mark[k];
1247 if (iMark) {
1248 // something in chunk - do all (as imark may change)
1249 int iLast = k<<CHECK_SHIFT;
1250 for (int i = iLast+BITS_PER_CHECK-1; i >= iLast; i-- ) {
1251 CoinFactorizationDouble pivotValue = region[i];
1252 if ( fabs ( pivotValue ) > tolerance ) {
1253 regionIndex[numberNonZero++] = i;
1254 CoinBigIndex j;
1255 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
1256 int iRow = column[j];
1257 CoinFactorizationDouble value = element[j];
1258 int iWord = iRow>>CHECK_SHIFT;
1259 int iBit = iRow-(iWord<<CHECK_SHIFT);
1260 if (mark[iWord]) {
1261 mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
1262 } else {
1263 mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
1264 }
1265 region[iRow] -= pivotValue*value;
1266 }
1267 } else {
1268 region[i] = 0.0;
1269 }
1270 }
1271 mark[k]=0;
1272 }
1273 }
1274#ifdef COIN_DEBUG
1275 for (int i=0;i<maximumRowsExtra_;i++) {
1276 assert (!mark[i]);
1277 }
1278#endif
1279 //set counts
1280 regionSparse->setNumElements ( numberNonZero );
1281}
1282/* updateColumnTransposeLSparse.
1283 Updates part of column transpose (BTRANL) sparse */
1284void
1285CoinFactorization::updateColumnTransposeLSparse
1286 ( CoinIndexedVector * regionSparse ) const
1287{
1288 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1289 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
1290 int numberNonZero = regionSparse->getNumElements ( );
1291 double tolerance = zeroTolerance_;
1292
1293 // use row copy of L
1294 const CoinFactorizationDouble * element = elementByRowL_.array();
1295 const CoinBigIndex * startRow = startRowL_.array();
1296 const int * column = indexColumnL_.array();
1297 // use sparse_ as temporary area
1298 // mark known to be zero
1299 int * COIN_RESTRICT stack = sparse_.array(); /* pivot */
1300 int * COIN_RESTRICT list = stack + maximumRowsExtra_; /* final list */
1301 CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_); /* jnext */
1302 char * COIN_RESTRICT mark = reinterpret_cast<char *> (next + maximumRowsExtra_);
1303 int nList;
1304 int number = numberNonZero;
1305#ifdef COIN_DEBUG
1306 for (i=0;i<maximumRowsExtra_;i++) {
1307 assert (!mark[i]);
1308 }
1309#endif
1310 nList=0;
1311 for (int k=0;k<number;k++) {
1312 int kPivot=regionIndex[k];
1313 if(!mark[kPivot]&&region[kPivot]) {
1314 stack[0]=kPivot;
1315 CoinBigIndex j=startRow[kPivot+1]-1;
1316 int nStack=0;
1317 while (nStack>=0) {
1318 /* take off stack */
1319 if (j>=startRow[kPivot]) {
1320 int jPivot=column[j--];
1321 /* put back on stack */
1322 next[nStack] =j;
1323 if (!mark[jPivot]) {
1324 /* and new one */
1325 kPivot=jPivot;
1326 j = startRow[kPivot+1]-1;
1327 stack[++nStack]=kPivot;
1328 mark[kPivot]=1;
1329 next[nStack]=j;
1330 }
1331 } else {
1332 /* finished so mark */
1333 list[nList++]=kPivot;
1334 mark[kPivot]=1;
1335 --nStack;
1336 if (nStack>=0) {
1337 kPivot=stack[nStack];
1338 j=next[nStack];
1339 }
1340 }
1341 }
1342 }
1343 }
1344 numberNonZero=0;
1345 for (int i=nList-1;i>=0;i--) {
1346 int iPivot = list[i];
1347 mark[iPivot]=0;
1348 CoinFactorizationDouble pivotValue = region[iPivot];
1349 if ( fabs ( pivotValue ) > tolerance ) {
1350 regionIndex[numberNonZero++] = iPivot;
1351 CoinBigIndex j;
1352 for ( j = startRow[iPivot]; j < startRow[iPivot+1]; j ++ ) {
1353 int iRow = column[j];
1354 CoinFactorizationDouble value = element[j];
1355 region[iRow] -= value * pivotValue;
1356 }
1357 } else {
1358 region[iPivot]=0.0;
1359 }
1360 }
1361 //set counts
1362 regionSparse->setNumElements ( numberNonZero );
1363}
1364// updateColumnTransposeL. Updates part of column transpose (BTRANL)
1365void
1366CoinFactorization::updateColumnTransposeL ( CoinIndexedVector * regionSparse ) const
1367{
1368 int number = regionSparse->getNumElements ( );
1369 if (!numberL_&&!numberDense_) {
1370 if (sparse_.array()||number<numberRows_)
1371 return;
1372 }
1373 int goSparse;
1374 // Guess at number at end
1375 // we may need to rethink on dense
1376 if (sparseThreshold_>0) {
1377 if (btranAverageAfterL_) {
1378 int newNumber = static_cast<int> (number*btranAverageAfterL_);
1379 if (newNumber< sparseThreshold_)
1380 goSparse = 2;
1381 else if (newNumber< sparseThreshold2_)
1382 goSparse = 1;
1383 else
1384 goSparse = 0;
1385 } else {
1386 if (number<sparseThreshold_)
1387 goSparse = 2;
1388 else
1389 goSparse = 0;
1390 }
1391 } else {
1392 goSparse=-1;
1393 }
1394#ifdef DENSE_CODE
1395 if (numberDense_) {
1396 //take off list
1397 int lastSparse = numberRows_-numberDense_;
1398 int number = regionSparse->getNumElements();
1399 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1400 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
1401 int i=0;
1402 bool doDense=false;
1403 if (number<=numberRows_) {
1404 while (i<number) {
1405 int iRow = regionIndex[i];
1406 if (iRow>=lastSparse) {
1407 doDense=true;
1408 regionIndex[i] = regionIndex[--number];
1409 } else {
1410 i++;
1411 }
1412 }
1413 } else {
1414 for (i=numberRows_-1;i>=lastSparse;i--) {
1415 if (region[i]) {
1416 doDense=true;
1417 // numbers are all wrong - do a scan
1418 regionSparse->setNumElements(0);
1419 regionSparse->scan(0,lastSparse,zeroTolerance_);
1420 number=regionSparse->getNumElements();
1421 break;
1422 }
1423 }
1424 if (sparseThreshold_)
1425 goSparse=0;
1426 else
1427 goSparse=-1;
1428 }
1429 if (doDense) {
1430 regionSparse->setNumElements(number);
1431 char trans = 'T';
1432 int ione=1;
1433 int info;
1434 F77_FUNC(dgetrs,DGETRS)(&trans,&numberDense_,&ione,denseArea_,&numberDense_,
1435 densePermute_,region+lastSparse,&numberDense_,&info,1);
1436 //and scan again
1437 if (goSparse>0||!numberL_)
1438 regionSparse->scan(lastSparse,numberRows_,zeroTolerance_);
1439 }
1440 if (!numberL_) {
1441 // could be odd combination of sparse/dense
1442 if (number>numberRows_) {
1443 regionSparse->setNumElements(0);
1444 regionSparse->scan(0,numberRows_,zeroTolerance_);
1445 }
1446 return;
1447 }
1448 }
1449#endif
1450 switch (goSparse) {
1451 case -1: // No row copy
1452 updateColumnTransposeLDensish(regionSparse);
1453 break;
1454 case 0: // densish but by row
1455 updateColumnTransposeLByRow(regionSparse);
1456 break;
1457 case 1: // middling(and by row)
1458 updateColumnTransposeLSparsish(regionSparse);
1459 break;
1460 case 2: // sparse
1461 updateColumnTransposeLSparse(regionSparse);
1462 break;
1463 }
1464}
1465#if COIN_ONE_ETA_COPY
1466/* Combines BtranU and delete elements
1467 If deleted is NULL then delete elements
1468 otherwise store where elements are
1469*/
1470void
1471CoinFactorization::replaceColumnU ( CoinIndexedVector * regionSparse,
1472 CoinBigIndex * deleted,
1473 int internalPivotRow)
1474{
1475 double * COIN_RESTRICT region = regionSparse->denseVector();
1476 int * COIN_RESTRICT regionIndex = regionSparse->getIndices();
1477 double tolerance = zeroTolerance_;
1478 const CoinBigIndex *startColumn = startColumnU_.array();
1479 const int *indexRow = indexRowU_.array();
1480 CoinFactorizationDouble *element = elementU_.array();
1481 int numberNonZero = 0;
1482 const int *numberInColumn = numberInColumn_.array();
1483 //const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
1484 bool deleteNow=true;
1485 if (deleted) {
1486 deleteNow = false;
1487 deleted ++;
1488 }
1489 int nPut=0;
1490 for (int i = CoinMax(numberSlacks_,internalPivotRow) ;
1491 i < numberU_; i++ ) {
1492 assert (!region[i]);
1493 CoinFactorizationDouble pivotValue = 0.0; //region[i]*pivotRegion[i];
1494 //printf("Epv %g reg %g pr %g\n",
1495 // pivotValue,region[i],pivotRegion[i]);
1496 //pivotValue = region[i];
1497 for (CoinBigIndex j= startColumn[i] ;
1498 j < startColumn[i]+numberInColumn[i]; j++ ) {
1499 int iRow = indexRow[j];
1500 CoinFactorizationDouble value = element[j];
1501 if (iRow!=internalPivotRow) {
1502 pivotValue -= value * region[iRow];
1503 } else {
1504 assert (!region[iRow]);
1505 pivotValue += value;
1506 if (deleteNow)
1507 element[j]=0.0;
1508 else
1509 deleted[nPut++]=j;
1510 }
1511 }
1512 if ( fabs ( pivotValue ) > tolerance ) {
1513 regionIndex[numberNonZero++] = i;
1514 region[i] = pivotValue;
1515 } else {
1516 region[i] = 0;
1517 }
1518 }
1519 if (!deleteNow) {
1520 deleted--;
1521 deleted[0]=nPut;
1522 }
1523 regionSparse->setNumElements(numberNonZero);
1524}
1525/* Updates part of column transpose (BTRANU) by column
1526 assumes index is sorted i.e. region is correct */
1527void
1528CoinFactorization::updateColumnTransposeUByColumn ( CoinIndexedVector * regionSparse,
1529 int smallestIndex) const
1530{
1531 //CoinIndexedVector temp = *regionSparse;
1532 //updateColumnTransposeUDensish(&temp,smallestIndex);
1533 double * COIN_RESTRICT region = regionSparse->denseVector();
1534 int * COIN_RESTRICT regionIndex = regionSparse->getIndices();
1535 double tolerance = zeroTolerance_;
1536 const CoinBigIndex *startColumn = startColumnU_.array();
1537 const int *indexRow = indexRowU_.array();
1538 const CoinFactorizationDouble *element = elementU_.array();
1539 int numberNonZero = 0;
1540 const int *numberInColumn = numberInColumn_.array();
1541 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
1542
1543 for (int i = smallestIndex;i<numberSlacks_; i++) {
1544 double value = region[i];
1545 if ( value ) {
1546 //region[i]=-value;
1547 regionIndex[numberNonZero]=i;
1548 if ( fabs(value) > tolerance )
1549 numberNonZero++;
1550 else
1551 region[i]=0.0;
1552 }
1553 }
1554 for (int i = CoinMax(numberSlacks_,smallestIndex) ;
1555 i < numberU_; i++ ) {
1556 CoinFactorizationDouble pivotValue = region[i]*pivotRegion[i];
1557 //printf("pv %g reg %g pr %g\n",
1558 // pivotValue,region[i],pivotRegion[i]);
1559 pivotValue = region[i];
1560 for (CoinBigIndex j= startColumn[i] ;
1561 j < startColumn[i]+numberInColumn[i]; j++ ) {
1562 int iRow = indexRow[j];
1563 CoinFactorizationDouble value = element[j];
1564 pivotValue -= value * region[iRow];
1565 }
1566 if ( fabs ( pivotValue ) > tolerance ) {
1567 regionIndex[numberNonZero++] = i;
1568 region[i] = pivotValue;
1569 } else {
1570 region[i] = 0;
1571 }
1572 }
1573 regionSparse->setNumElements(numberNonZero);
1574 //double * region2 = temp.denseVector();
1575 //for (i=0;i<maximumRowsExtra_;i++) {
1576 //assert(fabs(region[i]-region2[i])<1.0e-4);
1577 //}
1578}
1579#endif
1580// Updates part of column transpose (BTRANR) when dense
1581void
1582CoinFactorization::updateColumnTransposeRDensish
1583( CoinIndexedVector * regionSparse ) const
1584{
1585 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1586
1587#ifdef COIN_DEBUG
1588 regionSparse->checkClean();
1589#endif
1590 int last = numberRowsExtra_-1;
1591
1592
1593 const int *indexRow = indexRowR_;
1594 const CoinFactorizationDouble *element = elementR_;
1595 const CoinBigIndex * startColumn = startColumnR_.array()-numberRows_;
1596 //move using permute_ (stored in inverse fashion)
1597 const int * permute = permute_.array();
1598
1599 for (int i = last ; i >= numberRows_; i-- ) {
1600 int putRow = permute[i];
1601 CoinFactorizationDouble pivotValue = region[i];
1602 //zero out old permuted
1603 region[i] = 0.0;
1604 if ( pivotValue ) {
1605 for (CoinBigIndex j = startColumn[i]; j < startColumn[i+1]; j++ ) {
1606 CoinFactorizationDouble value = element[j];
1607 int iRow = indexRow[j];
1608 region[iRow] -= value * pivotValue;
1609 }
1610 region[putRow] = pivotValue;
1611 //putRow must have been zero before so put on list ??
1612 //but can't catch up so will have to do L from end
1613 //unless we update lookBtran in replaceColumn - yes
1614 }
1615 }
1616}
1617// Updates part of column transpose (BTRANR) when sparse
1618void
1619CoinFactorization::updateColumnTransposeRSparse
1620( CoinIndexedVector * regionSparse ) const
1621{
1622 double * COIN_RESTRICT region = regionSparse->denseVector ( );
1623 int * COIN_RESTRICT regionIndex = regionSparse->getIndices ( );
1624 int numberNonZero = regionSparse->getNumElements ( );
1625 double tolerance = zeroTolerance_;
1626
1627#ifdef COIN_DEBUG
1628 regionSparse->checkClean();
1629#endif
1630 int last = numberRowsExtra_-1;
1631
1632
1633 const int *indexRow = indexRowR_;
1634 const CoinFactorizationDouble *element = elementR_;
1635 const CoinBigIndex * startColumn = startColumnR_.array()-numberRows_;
1636 //move using permute_ (stored in inverse fashion)
1637 const int * permute = permute_.array();
1638
1639 // we can use sparse_ as temporary array
1640 int * COIN_RESTRICT spare = sparse_.array();
1641 for (int i=0;i<numberNonZero;i++) {
1642 spare[regionIndex[i]]=i;
1643 }
1644 // still need to do in correct order (for now)
1645 for (int i = last ; i >= numberRows_; i-- ) {
1646 int putRow = permute[i];
1647 assert (putRow<=i);
1648 CoinFactorizationDouble pivotValue = region[i];
1649 //zero out old permuted
1650 region[i] = 0.0;
1651 if ( pivotValue ) {
1652 for (CoinBigIndex j = startColumn[i]; j < startColumn[i+1]; j++ ) {
1653 CoinFactorizationDouble value = element[j];
1654 int iRow = indexRow[j];
1655 CoinFactorizationDouble oldValue = region[iRow];
1656 CoinFactorizationDouble newValue = oldValue - value * pivotValue;
1657 if (oldValue) {
1658 if (newValue)
1659 region[iRow]=newValue;
1660 else
1661 region[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
1662 } else if (fabs(newValue)>tolerance) {
1663 region[iRow] = newValue;
1664 spare[iRow]=numberNonZero;
1665 regionIndex[numberNonZero++]=iRow;
1666 }
1667 }
1668 region[putRow] = pivotValue;
1669 // modify list
1670 int position=spare[i];
1671 regionIndex[position]=putRow;
1672 spare[putRow]=position;
1673 }
1674 }
1675 regionSparse->setNumElements(numberNonZero);
1676}
1677
1678// updateColumnTransposeR. Updates part of column (FTRANR)
1679void
1680CoinFactorization::updateColumnTransposeR ( CoinIndexedVector * regionSparse ) const
1681{
1682 if (numberRowsExtra_==numberRows_)
1683 return;
1684 int numberNonZero = regionSparse->getNumElements ( );
1685
1686 if (numberNonZero) {
1687 if (numberNonZero < (sparseThreshold_<<2)||(!numberL_&&sparse_.array())) {
1688 updateColumnTransposeRSparse ( regionSparse );
1689 if (collectStatistics_)
1690 btranCountAfterR_ += regionSparse->getNumElements();
1691 } else {
1692 updateColumnTransposeRDensish ( regionSparse );
1693 // we have lost indices
1694 // make sure won't try and go sparse again
1695 if (collectStatistics_)
1696 btranCountAfterR_ += CoinMin((numberNonZero<<1),numberRows_);
1697 regionSparse->setNumElements (numberRows_+1);
1698 }
1699 }
1700}
1701// makes a row copy of L
1702void
1703CoinFactorization::goSparse ( )
1704{
1705 if (!sparseThreshold_) {
1706 if (numberRows_>300) {
1707 if (numberRows_<10000) {
1708 sparseThreshold_=CoinMin(numberRows_/6,500);
1709 //sparseThreshold2_=sparseThreshold_;
1710 } else {
1711 sparseThreshold_=1000;
1712 //sparseThreshold2_=sparseThreshold_;
1713 }
1714 sparseThreshold2_=numberRows_>>2;
1715 } else {
1716 sparseThreshold_=0;
1717 sparseThreshold2_=0;
1718 }
1719 } else {
1720 if (!sparseThreshold_&&numberRows_>400) {
1721 sparseThreshold_=CoinMin((numberRows_-300)/9,1000);
1722 }
1723 sparseThreshold2_=sparseThreshold_;
1724 }
1725 if (!sparseThreshold_)
1726 return;
1727 // allow for stack, list, next and char map of mark
1728 int nRowIndex = (maximumRowsExtra_+CoinSizeofAsInt(int)-1)/
1729 CoinSizeofAsInt(char);
1730 int nInBig = static_cast<int>(sizeof(CoinBigIndex)/sizeof(int));
1731 assert (nInBig>=1);
1732 sparse_.conditionalNew( (2+nInBig)*maximumRowsExtra_ + nRowIndex );
1733 // zero out mark
1734 memset(sparse_.array()+(2+nInBig)*maximumRowsExtra_,
1735 0,maximumRowsExtra_*sizeof(char));
1736 elementByRowL_.conditionalDelete();
1737 indexColumnL_.conditionalDelete();
1738 startRowL_.conditionalNew(numberRows_+1);
1739 if (lengthAreaL_) {
1740 elementByRowL_.conditionalNew(lengthAreaL_);
1741 indexColumnL_.conditionalNew(lengthAreaL_);
1742 }
1743 // counts
1744 CoinBigIndex * COIN_RESTRICT startRowL = startRowL_.array();
1745 CoinZeroN(startRowL,numberRows_);
1746 const CoinBigIndex * startColumnL = startColumnL_.array();
1747 CoinFactorizationDouble * COIN_RESTRICT elementL = elementL_.array();
1748 const int * indexRowL = indexRowL_.array();
1749 for (int i=baseL_;i<baseL_+numberL_;i++) {
1750 for (CoinBigIndex j=startColumnL[i];j<startColumnL[i+1];j++) {
1751 int iRow = indexRowL[j];
1752 startRowL[iRow]++;
1753 }
1754 }
1755 // convert count to lasts
1756 CoinBigIndex count=0;
1757 for (int i=0;i<numberRows_;i++) {
1758 int numberInRow=startRowL[i];
1759 count += numberInRow;
1760 startRowL[i]=count;
1761 }
1762 startRowL[numberRows_]=count;
1763 // now insert
1764 CoinFactorizationDouble * COIN_RESTRICT elementByRowL = elementByRowL_.array();
1765 int * COIN_RESTRICT indexColumnL = indexColumnL_.array();
1766 for (int i=baseL_+numberL_-1;i>=baseL_;i--) {
1767 for (CoinBigIndex j=startColumnL[i];j<startColumnL[i+1];j++) {
1768 int iRow = indexRowL[j];
1769 CoinBigIndex start = startRowL[iRow]-1;
1770 startRowL[iRow]=start;
1771 elementByRowL[start]=elementL[j];
1772 indexColumnL[start]=i;
1773 }
1774 }
1775}
1776
1777// set sparse threshold
1778void
1779CoinFactorization::sparseThreshold ( int value )
1780{
1781 if (value>0&&sparseThreshold_) {
1782 sparseThreshold_=value;
1783 sparseThreshold2_= sparseThreshold_;
1784 } else if (!value&&sparseThreshold_) {
1785 // delete copy
1786 sparseThreshold_=0;
1787 sparseThreshold2_= 0;
1788 elementByRowL_.conditionalDelete();
1789 startRowL_.conditionalDelete();
1790 indexColumnL_.conditionalDelete();
1791 sparse_.conditionalDelete();
1792 } else if (value>0&&!sparseThreshold_) {
1793 if (value>1)
1794 sparseThreshold_=value;
1795 else
1796 sparseThreshold_=0;
1797 sparseThreshold2_= sparseThreshold_;
1798 goSparse();
1799 }
1800}
1801void CoinFactorization::maximumPivots ( int value )
1802{
1803 if (value>0) {
1804 maximumPivots_=value;
1805 }
1806}
1807void CoinFactorization::messageLevel ( int value )
1808{
1809 if (value>0&&value<16) {
1810 messageLevel_=value;
1811 }
1812}
1813void CoinFactorization::pivotTolerance ( double value )
1814{
1815 if (value>0.0&&value<=1.0) {
1816 pivotTolerance_=value;
1817 }
1818}
1819void CoinFactorization::zeroTolerance ( double value )
1820{
1821 if (value>0.0&&value<1.0) {
1822 zeroTolerance_=value;
1823 }
1824}
1825#ifndef COIN_FAST_CODE
1826void CoinFactorization::slackValue ( double value )
1827{
1828 if (value>=0.0) {
1829 slackValue_=1.0;
1830 } else {
1831 slackValue_=-1.0;
1832 }
1833}
1834#endif
1835// Reset all sparsity etc statistics
1836void CoinFactorization::resetStatistics()
1837{
1838 collectStatistics_=false;
1839
1840 /// Below are all to collect
1841 ftranCountInput_=0.0;
1842 ftranCountAfterL_=0.0;
1843 ftranCountAfterR_=0.0;
1844 ftranCountAfterU_=0.0;
1845 btranCountInput_=0.0;
1846 btranCountAfterU_=0.0;
1847 btranCountAfterR_=0.0;
1848 btranCountAfterL_=0.0;
1849
1850 /// We can roll over factorizations
1851 numberFtranCounts_=0;
1852 numberBtranCounts_=0;
1853
1854 /// While these are averages collected over last
1855 ftranAverageAfterL_=0.0;
1856 ftranAverageAfterR_=0.0;
1857 ftranAverageAfterU_=0.0;
1858 btranAverageAfterU_=0.0;
1859 btranAverageAfterR_=0.0;
1860 btranAverageAfterL_=0.0;
1861}
1862/* getColumnSpaceIterate. Gets space for one extra U element in Column
1863 may have to do compression (returns true)
1864 also moves existing vector.
1865 Returns -1 if no memory or where element was put
1866 Used by replaceRow (turns off R version) */
1867CoinBigIndex
1868CoinFactorization::getColumnSpaceIterate ( int iColumn, double value,
1869 int iRow)
1870{
1871 if (numberInColumnPlus_.array()) {
1872 numberInColumnPlus_.conditionalDelete();
1873 }
1874 int * COIN_RESTRICT numberInRow = numberInRow_.array();
1875 int * COIN_RESTRICT numberInColumn = numberInColumn_.array();
1876 int * COIN_RESTRICT nextColumn = nextColumn_.array();
1877 int * COIN_RESTRICT lastColumn = lastColumn_.array();
1878 int number = numberInColumn[iColumn];
1879 int iNext=nextColumn[iColumn];
1880 CoinBigIndex * COIN_RESTRICT startColumnU = startColumnU_.array();
1881 CoinBigIndex * COIN_RESTRICT startRowU = startRowU_.array();
1882 CoinBigIndex space = startColumnU[iNext]-startColumnU[iColumn];
1883 CoinBigIndex put;
1884 CoinBigIndex * COIN_RESTRICT convertRowToColumnU = convertRowToColumnU_.array();
1885 int * COIN_RESTRICT indexColumnU = indexColumnU_.array();
1886 CoinFactorizationDouble * COIN_RESTRICT elementU = elementU_.array();
1887 int * COIN_RESTRICT indexRowU = indexRowU_.array();
1888 if ( space < number + 1 ) {
1889 //see if it can go in at end
1890 if (lengthAreaU_-startColumnU[maximumColumnsExtra_]<number+1) {
1891 //compression
1892 int jColumn = nextColumn[maximumColumnsExtra_];
1893 CoinBigIndex put = 0;
1894 while ( jColumn != maximumColumnsExtra_ ) {
1895 //move
1896 CoinBigIndex get;
1897 CoinBigIndex getEnd;
1898
1899 get = startColumnU[jColumn];
1900 getEnd = get + numberInColumn[jColumn];
1901 startColumnU[jColumn] = put;
1902 CoinBigIndex i;
1903 for ( i = get; i < getEnd; i++ ) {
1904 CoinFactorizationDouble value = elementU[i];
1905 if (value) {
1906 indexRowU[put] = indexRowU[i];
1907 elementU[put] = value;
1908 put++;
1909 } else {
1910 numberInColumn[jColumn]--;
1911 }
1912 }
1913 jColumn = nextColumn[jColumn];
1914 }
1915 numberCompressions_++;
1916 startColumnU[maximumColumnsExtra_]=put;
1917 //space for cross reference
1918 CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
1919 CoinBigIndex j = 0;
1920 CoinBigIndex *startRow = startRowU_.array();
1921
1922 int iRow;
1923 for ( iRow = 0; iRow < numberRowsExtra_; iRow++ ) {
1924 startRow[iRow] = j;
1925 j += numberInRow[iRow];
1926 }
1927 factorElements_=j;
1928
1929 CoinZeroN ( numberInRow, numberRowsExtra_ );
1930 int i;
1931 for ( i = 0; i < numberRowsExtra_; i++ ) {
1932 CoinBigIndex start = startColumnU[i];
1933 CoinBigIndex end = start + numberInColumn[i];
1934
1935 CoinBigIndex j;
1936 for ( j = start; j < end; j++ ) {
1937 int iRow = indexRowU[j];
1938 int iLook = numberInRow[iRow];
1939
1940 numberInRow[iRow] = iLook + 1;
1941 CoinBigIndex k = startRow[iRow] + iLook;
1942
1943 indexColumnU[k] = i;
1944 convertRowToColumn[k] = j;
1945 }
1946 }
1947 }
1948 // Still may not be room (as iColumn was still in)
1949 if (lengthAreaU_-startColumnU[maximumColumnsExtra_]>=number+1) {
1950 int next = nextColumn[iColumn];
1951 int last = lastColumn[iColumn];
1952
1953 //out
1954 nextColumn[last] = next;
1955 lastColumn[next] = last;
1956
1957 put = startColumnU[maximumColumnsExtra_];
1958 //in at end
1959 last = lastColumn[maximumColumnsExtra_];
1960 nextColumn[last] = iColumn;
1961 lastColumn[maximumColumnsExtra_] = iColumn;
1962 lastColumn[iColumn] = last;
1963 nextColumn[iColumn] = maximumColumnsExtra_;
1964
1965 //move
1966 CoinBigIndex get = startColumnU[iColumn];
1967 startColumnU[iColumn] = put;
1968 int i = 0;
1969 for (i=0 ; i < number; i ++ ) {
1970 CoinFactorizationDouble value = elementU[get];
1971 int iRow=indexRowU[get++];
1972 if (value) {
1973 elementU[put]= value;
1974 int n=numberInRow[iRow];
1975 CoinBigIndex start = startRowU[iRow];
1976 CoinBigIndex j;
1977 for (j=start;j<start+n;j++) {
1978 if (indexColumnU[j]==iColumn) {
1979 convertRowToColumnU[j]=put;
1980 break;
1981 }
1982 }
1983 assert (j<start+n);
1984 indexRowU[put++] = iRow;
1985 } else {
1986 assert (!numberInRow[iRow]);
1987 numberInColumn[iColumn]--;
1988 }
1989 }
1990 //insert
1991 int n=numberInRow[iRow];
1992 CoinBigIndex start = startRowU[iRow];
1993 CoinBigIndex j;
1994 for (j=start;j<start+n;j++) {
1995 if (indexColumnU[j]==iColumn) {
1996 convertRowToColumnU[j]=put;
1997 break;
1998 }
1999 }
2000 assert (j<start+n);
2001 elementU[put]=value;
2002 indexRowU[put]=iRow;
2003 numberInColumn[iColumn]++;
2004 //add 4 for luck
2005 startColumnU[maximumColumnsExtra_] = CoinMin(static_cast<CoinBigIndex> (put + 4) ,lengthAreaU_);
2006 } else {
2007 // no room
2008 put=-1;
2009 }
2010 } else {
2011 // just slot in
2012 put=startColumnU[iColumn]+numberInColumn[iColumn];
2013 int n=numberInRow[iRow];
2014 CoinBigIndex start = startRowU[iRow];
2015 CoinBigIndex j;
2016 for (j=start;j<start+n;j++) {
2017 if (indexColumnU[j]==iColumn) {
2018 convertRowToColumnU[j]=put;
2019 break;
2020 }
2021 }
2022 assert (j<start+n);
2023 elementU[put]=value;
2024 indexRowU[put]=iRow;
2025 numberInColumn[iColumn]++;
2026 }
2027 return put;
2028}
2029/* Replaces one Row in basis,
2030 At present assumes just a singleton on row is in basis
2031 returns 0=OK, 1=Probably OK, 2=singular, 3 no space */
2032int
2033CoinFactorization::replaceRow ( int whichRow, int iNumberInRow,
2034 const int indicesColumn[], const double elements[] )
2035{
2036 if (!iNumberInRow)
2037 return 0;
2038 int next = nextRow_.array()[whichRow];
2039 int * numberInRow = numberInRow_.array();
2040#ifndef NDEBUG
2041 const int * numberInColumn = numberInColumn_.array();
2042#endif
2043 int numberNow = numberInRow[whichRow];
2044 const CoinBigIndex * startRowU = startRowU_.array();
2045 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
2046 CoinBigIndex start = startRowU[whichRow];
2047 CoinFactorizationDouble * elementU = elementU_.array();
2048 CoinBigIndex *convertRowToColumnU = convertRowToColumnU_.array();
2049 if (numberNow&&numberNow<100) {
2050 int ind[100];
2051 CoinMemcpyN(indexColumnU_.array()+start,numberNow,ind);
2052 int i;
2053 for (i=0;i<iNumberInRow;i++) {
2054 int jColumn=indicesColumn[i];
2055 int k;
2056 for (k=0;k<numberNow;k++) {
2057 if (ind[k]==jColumn) {
2058 ind[k]=-1;
2059 break;
2060 }
2061 }
2062 if (k==numberNow) {
2063 printf("new column %d not in current\n",jColumn);
2064 } else {
2065 k=convertRowToColumnU[start+k];
2066 CoinFactorizationDouble oldValue = elementU[k];
2067 CoinFactorizationDouble newValue = elements[i]*pivotRegion[jColumn];
2068 if (fabs(oldValue-newValue)>1.0e-3)
2069 printf("column %d, old value %g new %g (el %g, piv %g)\n",jColumn,oldValue,
2070 newValue,elements[i],pivotRegion[jColumn]);
2071 }
2072 }
2073 for (i=0;i<numberNow;i++) {
2074 if (ind[i]>=0)
2075 printf("current column %d not in new\n",ind[i]);
2076 }
2077 assert (numberNow==iNumberInRow);
2078 }
2079 assert (numberInColumn[whichRow]==0);
2080 assert (pivotRegion[whichRow]==1.0);
2081 CoinBigIndex space;
2082 int returnCode=0;
2083
2084 space = startRowU[next] - (start+iNumberInRow);
2085 if ( space < 0 ) {
2086 if (!getRowSpaceIterate ( whichRow, iNumberInRow))
2087 returnCode=3;
2088 }
2089 //return 0;
2090 if (!returnCode) {
2091 int * indexColumnU = indexColumnU_.array();
2092 numberInRow[whichRow]=iNumberInRow;
2093 start = startRowU[whichRow];
2094 int i;
2095 for (i=0;i<iNumberInRow;i++) {
2096 int iColumn = indicesColumn[i];
2097 indexColumnU[start+i]=iColumn;
2098 assert (iColumn>whichRow);
2099 CoinFactorizationDouble value = elements[i]*pivotRegion[iColumn];
2100#if 0
2101 int k;
2102 bool found=false;
2103 for (k=startColumnU[iColumn];k<startColumnU[iColumn]+numberInColumn[iColumn];k++) {
2104 if (indexRowU[k]==whichRow) {
2105 assert (fabs(elementU[k]-value)<1.0e-3);
2106 found=true;
2107 break;
2108 }
2109 }
2110#if 0
2111 assert (found);
2112#else
2113 if (found) {
2114 int number = numberInColumn[iColumn]-1;
2115 numberInColumn[iColumn]=number;
2116 CoinBigIndex j=startColumnU[iColumn]+number;
2117 if (k<j) {
2118 int iRow=indexRowU[j];
2119 indexRowU[k]=iRow;
2120 elementU[k]=elementU[j];
2121 int n=numberInRow[iRow];
2122 CoinBigIndex start = startRowU[iRow];
2123 for (j=start;j<start+n;j++) {
2124 if (indexColumnU[j]==iColumn) {
2125 convertRowToColumnU[j]=k;
2126 break;
2127 }
2128 }
2129 assert (j<start+n);
2130 }
2131 }
2132 found=false;
2133#endif
2134 if (!found) {
2135#endif
2136 CoinBigIndex iWhere = getColumnSpaceIterate(iColumn,value,whichRow);
2137 if (iWhere>=0) {
2138 convertRowToColumnU[start+i] = iWhere;
2139 } else {
2140 returnCode=3;
2141 break;
2142 }
2143#if 0
2144 } else {
2145 convertRowToColumnU[start+i] = k;
2146 }
2147#endif
2148 }
2149 }
2150 return returnCode;
2151}
2152// Takes out all entries for given rows
2153void
2154CoinFactorization::emptyRows(int numberToEmpty, const int which[])
2155{
2156 int i;
2157 int * delRow = new int [maximumRowsExtra_];
2158 int * indexRowU = indexRowU_.array();
2159#ifndef NDEBUG
2160 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
2161#endif
2162 for (i=0;i<maximumRowsExtra_;i++)
2163 delRow[i]=0;
2164 int * numberInRow = numberInRow_.array();
2165 int * numberInColumn = numberInColumn_.array();
2166 CoinFactorizationDouble * elementU = elementU_.array();
2167 const CoinBigIndex * startColumnU = startColumnU_.array();
2168 for (i=0;i<numberToEmpty;i++) {
2169 int iRow = which[i];
2170 delRow[iRow]=1;
2171 assert (numberInColumn[iRow]==0);
2172 assert (pivotRegion[iRow]==1.0);
2173 numberInRow[iRow]=0;
2174 }
2175 for (i=0;i<numberU_;i++) {
2176 CoinBigIndex k;
2177 CoinBigIndex j=startColumnU[i];
2178 for (k=startColumnU[i];k<startColumnU[i]+numberInColumn[i];k++) {
2179 int iRow=indexRowU[k];
2180 if (!delRow[iRow]) {
2181 indexRowU[j]=indexRowU[k];
2182 elementU[j++]=elementU[k];
2183 }
2184 }
2185 numberInColumn[i] = j-startColumnU[i];
2186 }
2187 delete [] delRow;
2188 //space for cross reference
2189 CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
2190 CoinBigIndex j = 0;
2191 CoinBigIndex *startRow = startRowU_.array();
2192
2193 int iRow;
2194 for ( iRow = 0; iRow < numberRows_; iRow++ ) {
2195 startRow[iRow] = j;
2196 j += numberInRow[iRow];
2197 }
2198 factorElements_=j;
2199
2200 CoinZeroN ( numberInRow, numberRows_ );
2201
2202 int * indexColumnU = indexColumnU_.array();
2203 for ( i = 0; i < numberRows_; i++ ) {
2204 CoinBigIndex start = startColumnU[i];
2205 CoinBigIndex end = start + numberInColumn[i];
2206
2207 CoinBigIndex j;
2208 for ( j = start; j < end; j++ ) {
2209 int iRow = indexRowU[j];
2210 int iLook = numberInRow[iRow];
2211
2212 numberInRow[iRow] = iLook + 1;
2213 CoinBigIndex k = startRow[iRow] + iLook;
2214
2215 indexColumnU[k] = i;
2216 convertRowToColumn[k] = j;
2217 }
2218 }
2219}
2220// Updates part of column PFI (FTRAN)
2221void
2222CoinFactorization::updateColumnPFI ( CoinIndexedVector * regionSparse) const
2223{
2224 double *region = regionSparse->denseVector ( );
2225 int * regionIndex = regionSparse->getIndices();
2226 double tolerance = zeroTolerance_;
2227 const CoinBigIndex *startColumn = startColumnU_.array()+numberRows_;
2228 const int *indexRow = indexRowU_.array();
2229 const CoinFactorizationDouble *element = elementU_.array();
2230 int numberNonZero = regionSparse->getNumElements();
2231 int i;
2232#ifdef COIN_DEBUG
2233 for (i=0;i<numberNonZero;i++) {
2234 int iRow=regionIndex[i];
2235 assert (iRow>=0&&iRow<numberRows_);
2236 assert (region[iRow]);
2237 }
2238#endif
2239 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array()+numberRows_;
2240 const int *pivotColumn = pivotColumn_.array()+numberRows_;
2241
2242 for (i = 0 ; i <numberPivots_; i++ ) {
2243 int pivotRow=pivotColumn[i];
2244 CoinFactorizationDouble pivotValue = region[pivotRow];
2245 if (pivotValue) {
2246 if ( fabs ( pivotValue ) > tolerance ) {
2247 for (CoinBigIndex j= startColumn[i] ; j < startColumn[i+1]; j++ ) {
2248 int iRow = indexRow[j];
2249 CoinFactorizationDouble oldValue = region[iRow];
2250 CoinFactorizationDouble value = oldValue - pivotValue*element[j];
2251 if (!oldValue) {
2252 if (fabs(value)>tolerance) {
2253 region[iRow]=value;
2254 regionIndex[numberNonZero++]=iRow;
2255 }
2256 } else {
2257 if (fabs(value)>tolerance) {
2258 region[iRow]=value;
2259 } else {
2260 region[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2261 }
2262 }
2263 }
2264 pivotValue *= pivotRegion[i];
2265 region[pivotRow]=pivotValue;
2266 } else if (pivotValue) {
2267 region[pivotRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2268 }
2269 }
2270 }
2271 regionSparse->setNumElements ( numberNonZero );
2272}
2273// Updates part of column transpose PFI (BTRAN),
2274
2275void
2276CoinFactorization::updateColumnTransposePFI ( CoinIndexedVector * regionSparse) const
2277{
2278 double *region = regionSparse->denseVector ( );
2279 int numberNonZero = regionSparse->getNumElements();
2280 int *index = regionSparse->getIndices ( );
2281 int i;
2282#ifdef COIN_DEBUG
2283 for (i=0;i<numberNonZero;i++) {
2284 int iRow=index[i];
2285 assert (iRow>=0&&iRow<numberRows_);
2286 assert (region[iRow]);
2287 }
2288#endif
2289 const int * pivotColumn = pivotColumn_.array()+numberRows_;
2290 const CoinFactorizationDouble *pivotRegion = pivotRegion_.array()+numberRows_;
2291 double tolerance = zeroTolerance_;
2292
2293 const CoinBigIndex *startColumn = startColumnU_.array()+numberRows_;
2294 const int *indexRow = indexRowU_.array();
2295 const CoinFactorizationDouble *element = elementU_.array();
2296
2297 for (i=numberPivots_-1 ; i>=0; i-- ) {
2298 int pivotRow = pivotColumn[i];
2299 CoinFactorizationDouble pivotValue = region[pivotRow]*pivotRegion[i];
2300 for (CoinBigIndex j= startColumn[i] ; j < startColumn[i+1]; j++ ) {
2301 int iRow = indexRow[j];
2302 CoinFactorizationDouble value = element[j];
2303 pivotValue -= value * region[iRow];
2304 }
2305 //pivotValue *= pivotRegion[i];
2306 if ( fabs ( pivotValue ) > tolerance ) {
2307 if (!region[pivotRow])
2308 index[numberNonZero++] = pivotRow;
2309 region[pivotRow] = pivotValue;
2310 } else {
2311 if (region[pivotRow])
2312 region[pivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
2313 }
2314 }
2315 //set counts
2316 regionSparse->setNumElements ( numberNonZero );
2317}
2318/* Replaces one Column to basis for PFI
2319 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
2320 Not sure what checkBeforeModifying means for PFI.
2321*/
2322int
2323CoinFactorization::replaceColumnPFI ( CoinIndexedVector * regionSparse,
2324 int pivotRow,
2325 double alpha)
2326{
2327 CoinBigIndex *startColumn=startColumnU_.array()+numberRows_;
2328 int *indexRow=indexRowU_.array();
2329 CoinFactorizationDouble *element=elementU_.array();
2330 CoinFactorizationDouble * pivotRegion = pivotRegion_.array()+numberRows_;
2331 // This has incoming column
2332 const double *region = regionSparse->denseVector ( );
2333 const int * index = regionSparse->getIndices();
2334 int numberNonZero = regionSparse->getNumElements();
2335
2336 int iColumn = numberPivots_;
2337
2338 if (!iColumn)
2339 startColumn[0]=startColumn[maximumColumnsExtra_];
2340 CoinBigIndex start = startColumn[iColumn];
2341
2342 //return at once if too many iterations
2343 if ( numberPivots_ >= maximumPivots_ ) {
2344 return 5;
2345 }
2346 if ( lengthAreaU_ - ( start + numberNonZero ) < 0) {
2347 return 3;
2348 }
2349
2350 int i;
2351 if (numberPivots_) {
2352 if (fabs(alpha)<1.0e-5) {
2353 if (fabs(alpha)<1.0e-7)
2354 return 2;
2355 else
2356 return 1;
2357 }
2358 } else {
2359 if (fabs(alpha)<1.0e-8)
2360 return 2;
2361 }
2362 CoinFactorizationDouble pivotValue = 1.0/alpha;
2363 pivotRegion[iColumn]=pivotValue;
2364 double tolerance = zeroTolerance_;
2365 const int * pivotColumn = pivotColumn_.array();
2366 // Operations done before permute back
2367 if (regionSparse->packedMode()) {
2368 for ( i = 0; i < numberNonZero; i++ ) {
2369 int iRow = index[i];
2370 if (iRow!=pivotRow) {
2371 if ( fabs ( region[i] ) > tolerance ) {
2372 indexRow[start]=pivotColumn[iRow];
2373 element[start++]=region[i]*pivotValue;
2374 }
2375 }
2376 }
2377 } else {
2378 for ( i = 0; i < numberNonZero; i++ ) {
2379 int iRow = index[i];
2380 if (iRow!=pivotRow) {
2381 if ( fabs ( region[iRow] ) > tolerance ) {
2382 indexRow[start]=pivotColumn[iRow];
2383 element[start++]=region[iRow]*pivotValue;
2384 }
2385 }
2386 }
2387 }
2388 numberPivots_++;
2389 numberNonZero=start-startColumn[iColumn];
2390 startColumn[numberPivots_]=start;
2391 totalElements_ += numberNonZero;
2392 int * pivotColumn2 = pivotColumn_.array()+numberRows_;
2393 pivotColumn2[iColumn]=pivotColumn[pivotRow];
2394 return 0;
2395}
2396// =
2397CoinFactorization & CoinFactorization::operator = ( const CoinFactorization & other ) {
2398 if (this != &other) {
2399 gutsOfDestructor();
2400 gutsOfInitialize(3);
2401 persistenceFlag_=other.persistenceFlag_;
2402 gutsOfCopy(other);
2403 }
2404 return *this;
2405}
2406void CoinFactorization::gutsOfCopy(const CoinFactorization &other)
2407{
2408 elementU_.allocate(other.elementU_, other.lengthAreaU_ *CoinSizeofAsInt(CoinFactorizationDouble));
2409 indexRowU_.allocate(other.indexRowU_, other.lengthAreaU_*CoinSizeofAsInt(int) );
2410 elementL_.allocate(other.elementL_, other.lengthAreaL_*CoinSizeofAsInt(CoinFactorizationDouble) );
2411 indexRowL_.allocate(other.indexRowL_, other.lengthAreaL_*CoinSizeofAsInt(int) );
2412 startColumnL_.allocate(other.startColumnL_,(other.numberRows_ + 1)*CoinSizeofAsInt(CoinBigIndex) );
2413 int extraSpace;
2414 if (other.numberInColumnPlus_.array()) {
2415 extraSpace = other.maximumPivots_ + 1 + other.maximumColumnsExtra_ + 1;
2416 } else {
2417 extraSpace = other.maximumPivots_ + 1 ;
2418 }
2419 startColumnR_.allocate(other.startColumnR_,extraSpace*CoinSizeofAsInt(CoinBigIndex));
2420 pivotRegion_.allocate(other.pivotRegion_, (other.maximumRowsExtra_ + 1 )*CoinSizeofAsInt(CoinFactorizationDouble));
2421 permuteBack_.allocate(other.permuteBack_,(other.maximumRowsExtra_ + 1)*CoinSizeofAsInt(int));
2422 permute_.allocate(other.permute_,(other.maximumRowsExtra_ + 1)*CoinSizeofAsInt(int));
2423 pivotColumnBack_.allocate(other.pivotColumnBack_,(other.maximumRowsExtra_ + 1)*CoinSizeofAsInt(int));
2424 firstCount_.allocate(other.firstCount_,(other.maximumRowsExtra_ + 1)*CoinSizeofAsInt(int));
2425 startColumnU_.allocate(other.startColumnU_, (other.maximumColumnsExtra_ + 1 )*CoinSizeofAsInt(CoinBigIndex));
2426 numberInColumn_.allocate(other.numberInColumn_, (other.maximumColumnsExtra_ + 1 )*CoinSizeofAsInt(int));
2427 pivotColumn_.allocate(other.pivotColumn_,(other.maximumColumnsExtra_ + 1)*CoinSizeofAsInt(int));
2428 nextColumn_.allocate(other.nextColumn_, (other.maximumColumnsExtra_ + 1 )*CoinSizeofAsInt(int));
2429 lastColumn_.allocate(other.lastColumn_, (other.maximumColumnsExtra_ + 1 )*CoinSizeofAsInt(int));
2430 indexColumnU_.allocate(other.indexColumnU_, other.lengthAreaU_*CoinSizeofAsInt(int) );
2431 nextRow_.allocate(other.nextRow_,(other.maximumRowsExtra_ + 1)*CoinSizeofAsInt(int));
2432 lastRow_.allocate( other.lastRow_,(other.maximumRowsExtra_ + 1 )*CoinSizeofAsInt(int));
2433 const CoinBigIndex * convertUOther = other.convertRowToColumnU_.array();
2434#if COIN_ONE_ETA_COPY
2435 if (convertUOther) {
2436#endif
2437 convertRowToColumnU_.allocate(other.convertRowToColumnU_, other.lengthAreaU_*CoinSizeofAsInt(CoinBigIndex) );
2438 startRowU_.allocate(other.startRowU_,(other.maximumRowsExtra_ + 1)*CoinSizeofAsInt(CoinBigIndex));
2439 numberInRow_.allocate(other.numberInRow_, (other.maximumRowsExtra_ + 1 )*CoinSizeofAsInt(int));
2440#if COIN_ONE_ETA_COPY
2441 }
2442#endif
2443 if (other.sparseThreshold_) {
2444 elementByRowL_.allocate(other.elementByRowL_, other.lengthAreaL_ );
2445 indexColumnL_.allocate(other.indexColumnL_, other.lengthAreaL_ );
2446 startRowL_.allocate(other.startRowL_,other.numberRows_+1);
2447 }
2448 numberTrials_ = other.numberTrials_;
2449 biggerDimension_ = other.biggerDimension_;
2450 relaxCheck_ = other.relaxCheck_;
2451 numberSlacks_ = other.numberSlacks_;
2452 numberU_ = other.numberU_;
2453 maximumU_=other.maximumU_;
2454 lengthU_ = other.lengthU_;
2455 lengthAreaU_ = other.lengthAreaU_;
2456 numberL_ = other.numberL_;
2457 baseL_ = other.baseL_;
2458 lengthL_ = other.lengthL_;
2459 lengthAreaL_ = other.lengthAreaL_;
2460 numberR_ = other.numberR_;
2461 lengthR_ = other.lengthR_;
2462 lengthAreaR_ = other.lengthAreaR_;
2463 pivotTolerance_ = other.pivotTolerance_;
2464 zeroTolerance_ = other.zeroTolerance_;
2465#ifndef COIN_FAST_CODE
2466 slackValue_ = other.slackValue_;
2467#endif
2468 areaFactor_ = other.areaFactor_;
2469 numberRows_ = other.numberRows_;
2470 numberRowsExtra_ = other.numberRowsExtra_;
2471 maximumRowsExtra_ = other.maximumRowsExtra_;
2472 numberColumns_ = other.numberColumns_;
2473 numberColumnsExtra_ = other.numberColumnsExtra_;
2474 maximumColumnsExtra_ = other.maximumColumnsExtra_;
2475 maximumPivots_=other.maximumPivots_;
2476 numberGoodU_ = other.numberGoodU_;
2477 numberGoodL_ = other.numberGoodL_;
2478 numberPivots_ = other.numberPivots_;
2479 messageLevel_ = other.messageLevel_;
2480 totalElements_ = other.totalElements_;
2481 factorElements_ = other.factorElements_;
2482 status_ = other.status_;
2483 doForrestTomlin_ = other.doForrestTomlin_;
2484 collectStatistics_=other.collectStatistics_;
2485 ftranCountInput_=other.ftranCountInput_;
2486 ftranCountAfterL_=other.ftranCountAfterL_;
2487 ftranCountAfterR_=other.ftranCountAfterR_;
2488 ftranCountAfterU_=other.ftranCountAfterU_;
2489 btranCountInput_=other.btranCountInput_;
2490 btranCountAfterU_=other.btranCountAfterU_;
2491 btranCountAfterR_=other.btranCountAfterR_;
2492 btranCountAfterL_=other.btranCountAfterL_;
2493 numberFtranCounts_=other.numberFtranCounts_;
2494 numberBtranCounts_=other.numberBtranCounts_;
2495 ftranAverageAfterL_=other.ftranAverageAfterL_;
2496 ftranAverageAfterR_=other.ftranAverageAfterR_;
2497 ftranAverageAfterU_=other.ftranAverageAfterU_;
2498 btranAverageAfterU_=other.btranAverageAfterU_;
2499 btranAverageAfterR_=other.btranAverageAfterR_;
2500 btranAverageAfterL_=other.btranAverageAfterL_;
2501 biasLU_=other.biasLU_;
2502 sparseThreshold_=other.sparseThreshold_;
2503 sparseThreshold2_=other.sparseThreshold2_;
2504 CoinBigIndex space = lengthAreaL_ - lengthL_;
2505
2506 numberDense_ = other.numberDense_;
2507 denseThreshold_=other.denseThreshold_;
2508 if (numberDense_) {
2509 denseArea_ = new double [numberDense_*numberDense_];
2510 CoinMemcpyN(other.denseArea_,
2511 numberDense_*numberDense_,denseArea_);
2512 densePermute_ = new int [numberDense_];
2513 CoinMemcpyN(other.densePermute_,
2514 numberDense_,densePermute_);
2515 }
2516
2517 lengthAreaR_ = space;
2518 elementR_ = elementL_.array() + lengthL_;
2519 indexRowR_ = indexRowL_.array() + lengthL_;
2520 workArea_ = other.workArea_;
2521 workArea2_ = other.workArea2_;
2522 //now copy everything
2523 //assuming numberRowsExtra==numberColumnsExtra
2524 if (numberRowsExtra_) {
2525 if (convertUOther) {
2526 CoinMemcpyN ( other.startRowU_.array(), numberRowsExtra_ + 1, startRowU_.array() );
2527 CoinMemcpyN ( other.numberInRow_.array(), numberRowsExtra_ + 1, numberInRow_.array() );
2528 startRowU_.array()[maximumRowsExtra_] = other.startRowU_.array()[maximumRowsExtra_];
2529 }
2530 CoinMemcpyN ( other.pivotRegion_.array(), numberRowsExtra_ , pivotRegion_.array() );
2531 CoinMemcpyN ( other.permuteBack_.array(), numberRowsExtra_ + 1, permuteBack_.array() );
2532 CoinMemcpyN ( other.permute_.array(), numberRowsExtra_ + 1, permute_.array() );
2533 CoinMemcpyN ( other.pivotColumnBack_.array(), numberRowsExtra_ + 1, pivotColumnBack_.array() );
2534 CoinMemcpyN ( other.firstCount_.array(), numberRowsExtra_ + 1, firstCount_.array() );
2535 CoinMemcpyN ( other.startColumnU_.array(), numberRowsExtra_ + 1, startColumnU_.array() );
2536 CoinMemcpyN ( other.numberInColumn_.array(), numberRowsExtra_ + 1, numberInColumn_.array() );
2537 CoinMemcpyN ( other.pivotColumn_.array(), numberRowsExtra_ + 1, pivotColumn_.array() );
2538 CoinMemcpyN ( other.nextColumn_.array(), numberRowsExtra_ + 1, nextColumn_.array() );
2539 CoinMemcpyN ( other.lastColumn_.array(), numberRowsExtra_ + 1, lastColumn_.array() );
2540 CoinMemcpyN ( other.startColumnR_.array() , numberRowsExtra_ - numberColumns_ + 1,
2541 startColumnR_.array() );
2542 //extra one at end
2543 startColumnU_.array()[maximumColumnsExtra_] =
2544 other.startColumnU_.array()[maximumColumnsExtra_];
2545 nextColumn_.array()[maximumColumnsExtra_] = other.nextColumn_.array()[maximumColumnsExtra_];
2546 lastColumn_.array()[maximumColumnsExtra_] = other.lastColumn_.array()[maximumColumnsExtra_];
2547 CoinMemcpyN ( other.nextRow_.array(), numberRowsExtra_ + 1, nextRow_.array() );
2548 CoinMemcpyN ( other.lastRow_.array(), numberRowsExtra_ + 1, lastRow_.array() );
2549 nextRow_.array()[maximumRowsExtra_] = other.nextRow_.array()[maximumRowsExtra_];
2550 lastRow_.array()[maximumRowsExtra_] = other.lastRow_.array()[maximumRowsExtra_];
2551 }
2552 CoinMemcpyN ( other.elementR_, lengthR_, elementR_ );
2553 CoinMemcpyN ( other.indexRowR_, lengthR_, indexRowR_ );
2554 //row and column copies of U
2555 /* as elements of U may have been zeroed but column counts zero
2556 copy all elements */
2557 const CoinBigIndex * startColumnU = startColumnU_.array();
2558 const int * numberInColumn = numberInColumn_.array();
2559#ifndef NDEBUG
2560 int maxU=0;
2561 for (int iRow = 0; iRow < numberRowsExtra_; iRow++ ) {
2562 CoinBigIndex start = startColumnU[iRow];
2563 int numberIn = numberInColumn[iRow];
2564 maxU = CoinMax(maxU,start+numberIn);
2565 }
2566 assert (maximumU_>=maxU);
2567#endif
2568 CoinMemcpyN ( other.elementU_.array() , maximumU_, elementU_.array() );
2569#if COIN_ONE_ETA_COPY
2570 if (convertUOther) {
2571#endif
2572 const int * indexColumnUOther = other.indexColumnU_.array();
2573 CoinBigIndex * convertU = convertRowToColumnU_.array();
2574 int * indexColumnU = indexColumnU_.array();
2575 const CoinBigIndex * startRowU = startRowU_.array();
2576 const int * numberInRow = numberInRow_.array();
2577 for (int iRow = 0; iRow < numberRowsExtra_; iRow++ ) {
2578 //row
2579 CoinBigIndex start = startRowU[iRow];
2580 int numberIn = numberInRow[iRow];
2581
2582 CoinMemcpyN ( indexColumnUOther + start, numberIn, indexColumnU + start );
2583 CoinMemcpyN (convertUOther + start , numberIn, convertU + start );
2584 }
2585#if COIN_ONE_ETA_COPY
2586 }
2587#endif
2588 const int * indexRowUOther = other.indexRowU_.array();
2589 int * indexRowU = indexRowU_.array();
2590 for (int iRow = 0; iRow < numberRowsExtra_; iRow++ ) {
2591 //column
2592 CoinBigIndex start = startColumnU[iRow];
2593 int numberIn = numberInColumn[iRow];
2594 CoinMemcpyN ( indexRowUOther + start, numberIn, indexRowU + start );
2595 }
2596 // L is contiguous
2597 if (numberRows_)
2598 CoinMemcpyN ( other.startColumnL_.array(), numberRows_ + 1, startColumnL_.array() );
2599 CoinMemcpyN ( other.elementL_.array(), lengthL_, elementL_.array() );
2600 CoinMemcpyN ( other.indexRowL_.array(), lengthL_, indexRowL_.array() );
2601 if (other.sparseThreshold_) {
2602 goSparse();
2603 }
2604}
2605// See if worth going sparse
2606void
2607CoinFactorization::checkSparse()
2608{
2609 // See if worth going sparse and when
2610 if (numberFtranCounts_>100) {
2611 ftranCountInput_= CoinMax(ftranCountInput_,1.0);
2612 ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,1.0);
2613 ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,1.0);
2614 ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,1.0);
2615 if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
2616 btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,1.0);
2617 btranAverageAfterR_ = CoinMax(btranCountAfterR_/btranCountAfterU_,1.0);
2618 btranAverageAfterL_ = CoinMax(btranCountAfterL_/btranCountAfterR_,1.0);
2619 } else {
2620 // we have not done any useful btrans (values pass?)
2621 btranAverageAfterU_ = 1.0;
2622 btranAverageAfterR_ = 1.0;
2623 btranAverageAfterL_ = 1.0;
2624 }
2625 }
2626 // scale back
2627
2628 ftranCountInput_ *= 0.8;
2629 ftranCountAfterL_ *= 0.8;
2630 ftranCountAfterR_ *= 0.8;
2631 ftranCountAfterU_ *= 0.8;
2632 btranCountInput_ *= 0.8;
2633 btranCountAfterU_ *= 0.8;
2634 btranCountAfterR_ *= 0.8;
2635 btranCountAfterL_ *= 0.8;
2636}
2637// Condition number - product of pivots after factorization
2638double
2639CoinFactorization::conditionNumber() const
2640{
2641 double condition = 1.0;
2642 const CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
2643 for (int i=0;i<numberRows_;i++) {
2644 condition *= pivotRegion[i];
2645 }
2646 condition = CoinMax(fabs(condition),1.0e-50);
2647 return 1.0/condition;
2648}
2649#ifdef COIN_DEVELOP
2650extern double ncall_DZ;
2651extern double nrow_DZ;
2652extern double nslack_DZ;
2653extern double nU_DZ;
2654extern double nnz_DZ;
2655extern double nDone_DZ;
2656extern double ncall_SZ;
2657extern double nrow_SZ;
2658extern double nslack_SZ;
2659extern double nU_SZ;
2660extern double nnz_SZ;
2661extern double nDone_SZ;
2662void print_fac_stats()
2663{
2664 double mult = ncall_DZ ? 1.0/ncall_DZ : 1.0;
2665 printf("UDen called %g times, average rows %g, average slacks %g, average (U-S) %g average nnz in %g average ops %g\n",
2666 ncall_DZ,mult*nrow_DZ,mult*nslack_DZ,mult*(nU_DZ-nslack_DZ),mult*nnz_DZ,mult*nDone_DZ);
2667 ncall_DZ=0.0;
2668 nrow_DZ=0.0;
2669 nslack_DZ=0.0;
2670 nU_DZ=0.0;
2671 nnz_DZ=0.0;
2672 nDone_DZ=0.0;
2673 mult = ncall_SZ ? 1.0/ncall_SZ : 1.0;
2674 printf("USpars called %g times, average rows %g, average slacks %g, average (U-S) %g average nnz in %g average ops %g\n",
2675 ncall_SZ,mult*nrow_SZ,mult*nslack_SZ,mult*(nU_SZ-nslack_SZ),mult*nnz_SZ,mult*nDone_SZ);
2676 ncall_SZ=0.0;
2677 nrow_SZ=0.0;
2678 nslack_SZ=0.0;
2679 nU_SZ=0.0;
2680 nnz_SZ=0.0;
2681 nDone_SZ=0.0;
2682}
2683#endif
2684