1/* $Id: CoinFactorization1.cpp 1448 2011-06-19 15:34:41Z stefan $ */
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 "CoinFactorization.hpp"
15#include "CoinIndexedVector.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "CoinPackedMatrix.hpp"
18#include "CoinFinite.hpp"
19#include <stdio.h>
20//:class CoinFactorization. Deals with Factorization and Updates
21// CoinFactorization. Constructor
22CoinFactorization::CoinFactorization ( )
23{
24 persistenceFlag_=0;
25 gutsOfInitialize(7);
26}
27
28/// Copy constructor
29CoinFactorization::CoinFactorization ( const CoinFactorization &other)
30{
31 persistenceFlag_=0;
32 gutsOfInitialize(3);
33 persistenceFlag_=other.persistenceFlag_;
34 gutsOfCopy(other);
35}
36/// The real work of constructors etc
37void CoinFactorization::gutsOfDestructor(int )
38{
39 delete [] denseArea_;
40 delete [] densePermute_;
41
42 elementU_.conditionalDelete();
43 startRowU_.conditionalDelete();
44 convertRowToColumnU_.conditionalDelete();
45 indexRowU_.conditionalDelete();
46 indexColumnU_.conditionalDelete();
47 startColumnU_.conditionalDelete();
48 elementL_.conditionalDelete();
49 indexRowL_.conditionalDelete();
50 startColumnL_.conditionalDelete();
51 startColumnR_.conditionalDelete();
52 numberInRow_.conditionalDelete();
53 numberInColumn_.conditionalDelete();
54 numberInColumnPlus_.conditionalDelete();
55 pivotColumn_.conditionalDelete();
56 pivotColumnBack_.conditionalDelete();
57 firstCount_.conditionalDelete();
58 nextCount_.conditionalDelete();
59 lastCount_.conditionalDelete();
60 permute_.conditionalDelete();
61 permuteBack_.conditionalDelete();
62 nextColumn_.conditionalDelete();
63 lastColumn_.conditionalDelete();
64 nextRow_.conditionalDelete();
65 lastRow_.conditionalDelete();
66 saveColumn_.conditionalDelete();
67 markRow_.conditionalDelete();
68 pivotRowL_.conditionalDelete();
69 pivotRegion_.conditionalDelete();
70 elementByRowL_.conditionalDelete();
71 startRowL_.conditionalDelete();
72 indexColumnL_.conditionalDelete();
73 sparse_.conditionalDelete();
74 workArea_.conditionalDelete();
75 workArea2_.conditionalDelete();
76 numberCompressions_ = 0;
77 biggerDimension_ = 0;
78 numberRows_ = 0;
79 numberRowsExtra_ = 0;
80 maximumRowsExtra_ = 0;
81 numberColumns_ = 0;
82 numberColumnsExtra_ = 0;
83 maximumColumnsExtra_ = 0;
84 numberGoodU_ = 0;
85 numberGoodL_ = 0;
86 totalElements_ = 0;
87 factorElements_ = 0;
88 status_ = -1;
89 numberSlacks_ = 0;
90 numberU_ = 0;
91 maximumU_=0;
92 lengthU_ = 0;
93 lengthAreaU_ = 0;
94 numberL_ = 0;
95 baseL_ = 0;
96 lengthL_ = 0;
97 lengthAreaL_ = 0;
98 numberR_ = 0;
99 lengthR_ = 0;
100 lengthAreaR_ = 0;
101 denseArea_=NULL;
102 densePermute_=NULL;
103 // next two offsets but this makes cleaner
104 elementR_=NULL;
105 indexRowR_=NULL;
106 numberDense_=0;
107 //persistenceFlag_=0;
108 ////denseThreshold_=0;
109
110}
111// type - 1 bit tolerances etc, 2 rest
112void CoinFactorization::gutsOfInitialize(int type)
113{
114 if ((type&2)!=0) {
115 numberCompressions_ = 0;
116 biggerDimension_ = 0;
117 numberRows_ = 0;
118 numberRowsExtra_ = 0;
119 maximumRowsExtra_ = 0;
120 numberColumns_ = 0;
121 numberColumnsExtra_ = 0;
122 maximumColumnsExtra_ = 0;
123 numberGoodU_ = 0;
124 numberGoodL_ = 0;
125 totalElements_ = 0;
126 factorElements_ = 0;
127 status_ = -1;
128 numberPivots_ = 0;
129 numberSlacks_ = 0;
130 numberU_ = 0;
131 maximumU_=0;
132 lengthU_ = 0;
133 lengthAreaU_ = 0;
134 numberL_ = 0;
135 baseL_ = 0;
136 lengthL_ = 0;
137 lengthAreaL_ = 0;
138 numberR_ = 0;
139 lengthR_ = 0;
140 lengthAreaR_ = 0;
141 elementR_ = NULL;
142 indexRowR_ = NULL;
143 // always switch off sparse
144 sparseThreshold_=0;
145 sparseThreshold2_= 0;
146 denseArea_ = NULL;
147 densePermute_=NULL;
148 numberDense_=0;
149 if (!persistenceFlag_) {
150 workArea_=CoinFactorizationDoubleArrayWithLength();
151 workArea2_=CoinUnsignedIntArrayWithLength();
152 pivotColumn_=CoinIntArrayWithLength();
153 }
154 }
155 // after 2 because of persistenceFlag_
156 if ((type&1)!=0) {
157 areaFactor_ = 0.0;
158 pivotTolerance_ = 1.0e-1;
159 zeroTolerance_ = 1.0e-13;
160#ifndef COIN_FAST_CODE
161 slackValue_ = -1.0;
162#endif
163 messageLevel_=0;
164 maximumPivots_=200;
165 numberTrials_ = 4;
166 relaxCheck_=1.0;
167#if DENSE_CODE==1
168 denseThreshold_=31;
169 denseThreshold_=71;
170#else
171 denseThreshold_=0;
172#endif
173 biasLU_=2;
174 doForrestTomlin_=true;
175 persistenceFlag_=0;
176 }
177 if ((type&4)!=0) {
178 // we need to get 1 element arrays for any with length n+1 !!
179 startColumnL_.conditionalNew (1);
180 startColumnR_.conditionalNew( 1 );
181 startRowU_.conditionalNew(1);
182 numberInRow_.conditionalNew(1);
183 nextRow_.conditionalNew(1);
184 lastRow_.conditionalNew( 1 );
185 pivotRegion_.conditionalNew(1);
186 permuteBack_.conditionalNew(1);
187 permute_.conditionalNew(1);
188 pivotColumnBack_.conditionalNew(1);
189 startColumnU_.conditionalNew( 1 );
190 numberInColumn_.conditionalNew(1);
191 numberInColumnPlus_.conditionalNew(1);
192 pivotColumn_.conditionalNew(1);
193 nextColumn_.conditionalNew(1);
194 lastColumn_.conditionalNew(1);
195 collectStatistics_=false;
196
197 // Below are all to collect
198 ftranCountInput_=0.0;
199 ftranCountAfterL_=0.0;
200 ftranCountAfterR_=0.0;
201 ftranCountAfterU_=0.0;
202 btranCountInput_=0.0;
203 btranCountAfterU_=0.0;
204 btranCountAfterR_=0.0;
205 btranCountAfterL_=0.0;
206
207 // We can roll over factorizations
208 numberFtranCounts_=0;
209 numberBtranCounts_=0;
210
211 // While these are averages collected over last
212 ftranAverageAfterL_=0;
213 ftranAverageAfterR_=0;
214 ftranAverageAfterU_=0;
215 btranAverageAfterU_=0;
216 btranAverageAfterR_=0;
217 btranAverageAfterL_=0;
218#ifdef ZEROFAULT
219 startColumnL_.array()[0] = 0;
220 startColumnR_.array()[0] = 0;
221 startRowU_.array()[0] = 0;
222 numberInRow_.array()[0] = 0;
223 nextRow_.array()[0] = 0;
224 lastRow_.array()[0] = 0;
225 pivotRegion_.array()[0] = 0.0;
226 permuteBack_.array()[0] = 0;
227 permute_.array()[0] = 0;
228 pivotColumnBack_.array()[0] = 0;
229 startColumnU_.array()[0] = 0;
230 numberInColumn_.array()[0] = 0;
231 numberInColumnPlus_.array()[0] = 0;
232 pivotColumn_.array()[0] = 0;
233 nextColumn_.array()[0] = 0;
234 lastColumn_.array()[0] = 0;
235#endif
236 }
237}
238//Part of LP
239int CoinFactorization::factorize (
240 const CoinPackedMatrix & matrix,
241 int rowIsBasic[],
242 int columnIsBasic[],
243 double areaFactor )
244{
245 // maybe for speed will be better to leave as many regions as possible
246 gutsOfDestructor();
247 gutsOfInitialize(2);
248 // ? is this correct
249 //if (biasLU_==2)
250 //biasLU_=3;
251 if (areaFactor)
252 areaFactor_ = areaFactor;
253 const int * row = matrix.getIndices();
254 const CoinBigIndex * columnStart = matrix.getVectorStarts();
255 const int * columnLength = matrix.getVectorLengths();
256 const double * element = matrix.getElements();
257 int numberRows=matrix.getNumRows();
258 int numberColumns=matrix.getNumCols();
259 int numberBasic = 0;
260 CoinBigIndex numberElements=0;
261 int numberRowBasic=0;
262
263 // compute how much in basis
264
265 int i;
266
267 for (i=0;i<numberRows;i++) {
268 if (rowIsBasic[i]>=0)
269 numberRowBasic++;
270 }
271
272 numberBasic = numberRowBasic;
273
274 for (i=0;i<numberColumns;i++) {
275 if (columnIsBasic[i]>=0) {
276 numberBasic++;
277 numberElements += columnLength[i];
278 }
279 }
280 if ( numberBasic > numberRows ) {
281 return -2; // say too many in basis
282 }
283 numberElements = 3 * numberBasic + 3 * numberElements + 20000;
284 getAreas ( numberRows, numberBasic, numberElements,
285 2 * numberElements );
286 //fill
287 //copy
288 numberBasic=0;
289 numberElements=0;
290 int * indexColumnU = indexColumnU_.array();
291 int * indexRowU = indexRowU_.array();
292 CoinFactorizationDouble * elementU = elementU_.array();
293 for (i=0;i<numberRows;i++) {
294 if (rowIsBasic[i]>=0) {
295 indexRowU[numberElements]=i;
296 indexColumnU[numberElements]=numberBasic;
297 elementU[numberElements++]=slackValue_;
298 numberBasic++;
299 }
300 }
301 for (i=0;i<numberColumns;i++) {
302 if (columnIsBasic[i]>=0) {
303 CoinBigIndex j;
304 for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
305 indexRowU[numberElements]=row[j];
306 indexColumnU[numberElements]=numberBasic;
307 elementU[numberElements++]=element[j];
308 }
309 numberBasic++;
310 }
311 }
312 lengthU_ = numberElements;
313 maximumU_ = numberElements;
314
315 preProcess ( 0 );
316 factor ( );
317 numberBasic=0;
318 if (status_ == 0) {
319 int * permuteBack = permuteBack_.array();
320 int * back = pivotColumnBack();
321 for (i=0;i<numberRows;i++) {
322 if (rowIsBasic[i]>=0) {
323 rowIsBasic[i]=permuteBack[back[numberBasic++]];
324 }
325 }
326 for (i=0;i<numberColumns;i++) {
327 if (columnIsBasic[i]>=0) {
328 columnIsBasic[i]=permuteBack[back[numberBasic++]];
329 }
330 }
331 // Set up permutation vector
332 // these arrays start off as copies of permute
333 // (and we could use permute_ instead of pivotColumn (not back though))
334 CoinMemcpyN ( permute_.array(), numberRows_ , pivotColumn_.array() );
335 CoinMemcpyN ( permuteBack_.array(), numberRows_ , pivotColumnBack() );
336 } else if (status_ == -1) {
337 const int * pivotColumn = pivotColumn_.array();
338 // mark as basic or non basic
339 for (i=0;i<numberRows_;i++) {
340 if (rowIsBasic[i]>=0) {
341 if (pivotColumn[numberBasic]>=0)
342 rowIsBasic[i]=pivotColumn[numberBasic];
343 else
344 rowIsBasic[i]=-1;
345 numberBasic++;
346 }
347 }
348 for (i=0;i<numberColumns;i++) {
349 if (columnIsBasic[i]>=0) {
350 if (pivotColumn[numberBasic]>=0)
351 columnIsBasic[i]=pivotColumn[numberBasic];
352 else
353 columnIsBasic[i]=-1;
354 numberBasic++;
355 }
356 }
357 }
358
359 return status_;
360}
361//Given as triplets
362int CoinFactorization::factorize (
363 int numberOfRows,
364 int numberOfColumns,
365 CoinBigIndex numberOfElements,
366 CoinBigIndex maximumL,
367 CoinBigIndex maximumU,
368 const int indicesRow[],
369 const int indicesColumn[],
370 const double elements[] ,
371 int permutation[],
372 double areaFactor)
373{
374 gutsOfDestructor();
375 gutsOfInitialize(2);
376 if (areaFactor)
377 areaFactor_ = areaFactor;
378 getAreas ( numberOfRows, numberOfColumns, maximumL, maximumU );
379 //copy
380 CoinMemcpyN ( indicesRow, numberOfElements, indexRowU_.array() );
381 CoinMemcpyN ( indicesColumn, numberOfElements, indexColumnU_.array() );
382 int i;
383 CoinFactorizationDouble * elementU = elementU_.array();
384 for (i=0;i<numberOfElements;i++)
385 elementU[i] = elements[i];
386 lengthU_ = numberOfElements;
387 maximumU_ = numberOfElements;
388 preProcess ( 0 );
389 factor ( );
390 //say which column is pivoting on which row
391 if (status_ == 0) {
392 int * permuteBack = permuteBack_.array();
393 int * back = pivotColumnBack();
394 // permute so slacks on own rows etc
395 for (i=0;i<numberOfColumns;i++) {
396 permutation[i]=permuteBack[back[i]];
397 }
398 // Set up permutation vector
399 // these arrays start off as copies of permute
400 // (and we could use permute_ instead of pivotColumn (not back though))
401 CoinMemcpyN ( permute_.array(), numberRows_ , pivotColumn_.array() );
402 CoinMemcpyN ( permuteBack_.array(), numberRows_ , pivotColumnBack() );
403 } else if (status_ == -1) {
404 const int * pivotColumn = pivotColumn_.array();
405 // mark as basic or non basic
406 for (i=0;i<numberOfColumns;i++) {
407 if (pivotColumn[i]>=0) {
408 permutation[i]=pivotColumn[i];
409 } else {
410 permutation[i]=-1;
411 }
412 }
413 }
414
415 return status_;
416}
417/* Two part version for flexibility
418 This part creates arrays for user to fill.
419 maximumL is guessed maximum size of L part of
420 final factorization, maximumU of U part. These are multiplied by
421 areaFactor which can be computed by user or internally.
422 returns 0 -okay, -99 memory */
423int
424CoinFactorization::factorizePart1 ( int numberOfRows,
425 int ,
426 CoinBigIndex numberOfElements,
427 int * indicesRow[],
428 int * indicesColumn[],
429 CoinFactorizationDouble * elements[],
430 double areaFactor)
431{
432 // maybe for speed will be better to leave as many regions as possible
433 gutsOfDestructor();
434 gutsOfInitialize(2);
435 if (areaFactor)
436 areaFactor_ = areaFactor;
437 CoinBigIndex numberElements = 3 * numberOfRows + 3 * numberOfElements + 20000;
438 getAreas ( numberOfRows, numberOfRows, numberElements,
439 2 * numberElements );
440 // need to trap memory for -99 code
441 *indicesRow = indexRowU_.array() ;
442 *indicesColumn = indexColumnU_.array() ;
443 *elements = elementU_.array() ;
444 lengthU_ = numberOfElements;
445 maximumU_ = numberElements;
446 return 0;
447}
448/* This is part two of factorization
449 Arrays belong to factorization and were returned by part 1
450 If status okay, permutation has pivot rows.
451 If status is singular, then basic variables have +1 and ones thrown out have -COIN_INT_MAX
452 to say thrown out.
453 returns 0 -okay, -1 singular, -99 memory */
454int
455CoinFactorization::factorizePart2 (int permutation[],int exactNumberElements)
456{
457 lengthU_ = exactNumberElements;
458 preProcess ( 0 );
459 factor ( );
460 //say which column is pivoting on which row
461 int i;
462 int * permuteBack = permuteBack_.array();
463 int * back = pivotColumnBack();
464 // permute so slacks on own rows etc
465 for (i=0;i<numberColumns_;i++) {
466 permutation[i]=permuteBack[back[i]];
467 }
468 if (status_ == 0) {
469 // Set up permutation vector
470 // these arrays start off as copies of permute
471 // (and we could use permute_ instead of pivotColumn (not back though))
472 CoinMemcpyN ( permute_.array(), numberRows_ , pivotColumn_.array() );
473 CoinMemcpyN ( permuteBack_.array(), numberRows_ , pivotColumnBack() );
474 } else if (status_ == -1) {
475 const int * pivotColumn = pivotColumn_.array();
476 // mark as basic or non basic
477 for (i=0;i<numberColumns_;i++) {
478 if (pivotColumn[i]>=0) {
479 permutation[i]=pivotColumn[i];
480 } else {
481 permutation[i]=-1;
482 }
483 }
484 }
485
486 return status_;
487}
488
489// ~CoinFactorization. Destructor
490CoinFactorization::~CoinFactorization ( )
491{
492 gutsOfDestructor();
493}
494
495// show_self. Debug show object
496void
497CoinFactorization::show_self ( ) const
498{
499 int i;
500
501 const int * pivotColumn = pivotColumn_.array();
502 for ( i = 0; i < numberRows_; i++ ) {
503 std::cout << "r " << i << " " << pivotColumn[i];
504 if (pivotColumnBack()) std::cout<< " " << pivotColumnBack()[i];
505 std::cout<< " " << permute_.array()[i];
506 if (permuteBack_.array()) std::cout<< " " << permuteBack_.array()[i];
507 std::cout<< " " << pivotRegion_.array()[i];
508 std::cout << std::endl;
509 }
510 for ( i = 0; i < numberRows_; i++ ) {
511 std::cout << "u " << i << " " << numberInColumn_.array()[i] << std::endl;
512 int j;
513 CoinSort_2(indexRowU_.array()+startColumnU_.array()[i],
514 indexRowU_.array()+startColumnU_.array()[i]+numberInColumn_.array()[i],
515 elementU_.array()+startColumnU_.array()[i]);
516 for ( j = startColumnU_.array()[i]; j < startColumnU_.array()[i] + numberInColumn_.array()[i];
517 j++ ) {
518 assert (indexRowU_.array()[j]>=0&&indexRowU_.array()[j]<numberRows_);
519 assert (elementU_.array()[j]>-1.0e50&&elementU_.array()[j]<1.0e50);
520 std::cout << indexRowU_.array()[j] << " " << elementU_.array()[j] << std::endl;
521 }
522 }
523 for ( i = 0; i < numberRows_; i++ ) {
524 std::cout << "l " << i << " " << startColumnL_.array()[i + 1] -
525 startColumnL_.array()[i] << std::endl;
526 CoinSort_2(indexRowL_.array()+startColumnL_.array()[i],
527 indexRowL_.array()+startColumnL_.array()[i+1],
528 elementL_.array()+startColumnL_.array()[i]);
529 int j;
530 for ( j = startColumnL_.array()[i]; j < startColumnL_.array()[i + 1]; j++ ) {
531 std::cout << indexRowL_.array()[j] << " " << elementL_.array()[j] << std::endl;
532 }
533 }
534
535}
536// sort so can compare
537void
538CoinFactorization::sort ( ) const
539{
540 int i;
541
542 for ( i = 0; i < numberRows_; i++ ) {
543 CoinSort_2(indexRowU_.array()+startColumnU_.array()[i],
544 indexRowU_.array()+startColumnU_.array()[i]+numberInColumn_.array()[i],
545 elementU_.array()+startColumnU_.array()[i]);
546 }
547 for ( i = 0; i < numberRows_; i++ ) {
548 CoinSort_2(indexRowL_.array()+startColumnL_.array()[i],
549 indexRowL_.array()+startColumnL_.array()[i+1],
550 elementL_.array()+startColumnL_.array()[i]);
551 }
552
553}
554
555// getAreas. Gets space for a factorization
556//called by constructors
557void
558CoinFactorization::getAreas ( int numberOfRows,
559 int numberOfColumns,
560 CoinBigIndex maximumL,
561 CoinBigIndex maximumU )
562{
563
564 numberRows_ = numberOfRows;
565 numberColumns_ = numberOfColumns;
566 maximumRowsExtra_ = numberRows_ + maximumPivots_;
567 numberRowsExtra_ = numberRows_;
568 maximumColumnsExtra_ = numberColumns_ + maximumPivots_;
569 numberColumnsExtra_ = numberColumns_;
570 lengthAreaU_ = maximumU;
571 lengthAreaL_ = maximumL;
572 if ( !areaFactor_ ) {
573 areaFactor_ = 1.0;
574 }
575 if ( areaFactor_ != 1.0 ) {
576 if ((messageLevel_&16)!=0)
577 printf("Increasing factorization areas by %g\n",areaFactor_);
578 lengthAreaU_ = static_cast<CoinBigIndex> (areaFactor_*lengthAreaU_);
579 lengthAreaL_ = static_cast<CoinBigIndex> (areaFactor_*lengthAreaL_);
580 }
581 elementU_.conditionalNew( lengthAreaU_ );
582 indexRowU_.conditionalNew( lengthAreaU_ );
583 indexColumnU_.conditionalNew( lengthAreaU_ );
584 elementL_.conditionalNew( lengthAreaL_ );
585 indexRowL_.conditionalNew( lengthAreaL_ );
586 if (persistenceFlag_) {
587 // But we can use all we have if bigger
588 int length;
589 length = CoinMin(elementU_.getSize(),indexRowU_.getSize());
590 if (length>lengthAreaU_) {
591 lengthAreaU_=length;
592 assert (indexColumnU_.getSize()==indexRowU_.getSize());
593 }
594 length = CoinMin(elementL_.getSize(),indexRowL_.getSize());
595 if (length>lengthAreaL_) {
596 lengthAreaL_=length;
597 }
598 }
599 startColumnL_.conditionalNew( numberRows_ + 1 );
600 startColumnL_.array()[0] = 0;
601 startRowU_.conditionalNew( maximumRowsExtra_ + 1);
602 // make sure this is valid
603 startRowU_.array()[maximumRowsExtra_]=0;
604 numberInRow_.conditionalNew( maximumRowsExtra_ + 1 );
605 markRow_.conditionalNew( numberRows_ );
606 pivotRowL_.conditionalNew( numberRows_ + 1 );
607 nextRow_.conditionalNew( maximumRowsExtra_ + 1 );
608 lastRow_.conditionalNew( maximumRowsExtra_ + 1 );
609 permute_.conditionalNew( maximumRowsExtra_ + 1 );
610 pivotRegion_.conditionalNew( maximumRowsExtra_ + 1 );
611#ifdef ZEROFAULT
612 memset(elementU_.array(),'a',lengthAreaU_*sizeof(CoinFactorizationDouble));
613 memset(indexRowU_.array(),'b',lengthAreaU_*sizeof(int));
614 memset(indexColumnU_.array(),'c',lengthAreaU_*sizeof(int));
615 memset(elementL_.array(),'d',lengthAreaL_*sizeof(CoinFactorizationDouble));
616 memset(indexRowL_.array(),'e',lengthAreaL_*sizeof(int));
617 memset(startColumnL_.array()+1,'f',numberRows_*sizeof(CoinBigIndex));
618 memset(startRowU_.array(),'g',maximumRowsExtra_*sizeof(CoinBigIndex));
619 memset(numberInRow_.array(),'h',(maximumRowsExtra_+1)*sizeof(int));
620 memset(markRow_.array(),'i',numberRows_*sizeof(int));
621 memset(pivotRowL_.array(),'j',(numberRows_+1)*sizeof(int));
622 memset(nextRow_.array(),'k',(maximumRowsExtra_+1)*sizeof(int));
623 memset(lastRow_.array(),'l',(maximumRowsExtra_+1)*sizeof(int));
624 memset(permute_.array(),'l',(maximumRowsExtra_+1)*sizeof(int));
625 memset(pivotRegion_.array(),'m',(maximumRowsExtra_+1)*sizeof(CoinFactorizationDouble));
626#endif
627 startColumnU_.conditionalNew( maximumColumnsExtra_ + 1 );
628 numberInColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
629 numberInColumnPlus_.conditionalNew( maximumColumnsExtra_ + 1 );
630 pivotColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
631 nextColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
632 lastColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
633 saveColumn_.conditionalNew( numberColumns_);
634#ifdef ZEROFAULT
635 memset(startColumnU_.array(),'a',(maximumColumnsExtra_+1)*sizeof(CoinBigIndex));
636 memset(numberInColumn_.array(),'b',(maximumColumnsExtra_+1)*sizeof(int));
637 memset(numberInColumnPlus_.array(),'c',(maximumColumnsExtra_+1)*sizeof(int));
638 memset(pivotColumn_.array(),'d',(maximumColumnsExtra_+1)*sizeof(int));
639 memset(nextColumn_.array(),'e',(maximumColumnsExtra_+1)*sizeof(int));
640 memset(lastColumn_.array(),'f',(maximumColumnsExtra_+1)*sizeof(int));
641#endif
642 if ( numberRows_ + numberColumns_ ) {
643 if ( numberRows_ > numberColumns_ ) {
644 biggerDimension_ = numberRows_;
645 } else {
646 biggerDimension_ = numberColumns_;
647 }
648 firstCount_.conditionalNew( CoinMax(biggerDimension_ + 2, maximumRowsExtra_+1) );
649 nextCount_.conditionalNew( numberRows_ + numberColumns_ );
650 lastCount_.conditionalNew( numberRows_ + numberColumns_ );
651#ifdef ZEROFAULT
652 memset(firstCount_.array(),'g',(biggerDimension_ + 2 )*sizeof(int));
653 memset(nextCount_.array(),'h',(numberRows_+numberColumns_)*sizeof(int));
654 memset(lastCount_.array(),'i',(numberRows_+numberColumns_)*sizeof(int));
655#endif
656 } else {
657 firstCount_.conditionalNew( 2 );
658 nextCount_.conditionalNew( 0 );
659 lastCount_.conditionalNew( 0 );
660#ifdef ZEROFAULT
661 memset(firstCount_.array(),'g', 2 *sizeof(int));
662#endif
663 biggerDimension_ = 0;
664 }
665}
666
667// preProcess. PreProcesses raw triplet data
668//state is 0 - triplets, 1 - some counts etc , 2 - ..
669void
670CoinFactorization::preProcess ( int state,
671 int )
672{
673 int *indexRow = indexRowU_.array();
674 int *indexColumn = indexColumnU_.array();
675 CoinFactorizationDouble *element = elementU_.array();
676 CoinBigIndex numberElements = lengthU_;
677 int *numberInRow = numberInRow_.array();
678 int *numberInColumn = numberInColumn_.array();
679 int *numberInColumnPlus = numberInColumnPlus_.array();
680 CoinBigIndex *startRow = startRowU_.array();
681 CoinBigIndex *startColumn = startColumnU_.array();
682 int numberRows = numberRows_;
683 int numberColumns = numberColumns_;
684
685 if (state<4)
686 totalElements_ = numberElements;
687 //state falls through to next state
688 switch ( state ) {
689 case 0: //counts
690 {
691 CoinZeroN ( numberInRow, numberRows + 1 );
692 CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 );
693 CoinBigIndex i;
694 for ( i = 0; i < numberElements; i++ ) {
695 int iRow = indexRow[i];
696 int iColumn = indexColumn[i];
697
698 numberInRow[iRow]++;
699 numberInColumn[iColumn]++;
700 }
701 }
702 case -1: //sort
703 case 1: //sort
704 {
705 CoinBigIndex i, k;
706
707 i = 0;
708 int iColumn;
709 for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
710 //position after end of Column
711 i += numberInColumn[iColumn];
712 startColumn[iColumn] = i;
713 }
714 for ( k = numberElements - 1; k >= 0; k-- ) {
715 int iColumn = indexColumn[k];
716
717 if ( iColumn >= 0 ) {
718 CoinFactorizationDouble value = element[k];
719 int iRow = indexRow[k];
720
721 indexColumn[k] = -1;
722 while ( true ) {
723 CoinBigIndex iLook = startColumn[iColumn] - 1;
724
725 startColumn[iColumn] = iLook;
726 CoinFactorizationDouble valueSave = element[iLook];
727 int iColumnSave = indexColumn[iLook];
728 int iRowSave = indexRow[iLook];
729
730 element[iLook] = value;
731 indexRow[iLook] = iRow;
732 indexColumn[iLook] = -1;
733 if ( iColumnSave >= 0 ) {
734 iColumn = iColumnSave;
735 value = valueSave;
736 iRow = iRowSave;
737 } else {
738 break;
739 }
740 } /* endwhile */
741 }
742 }
743 }
744 case 2: //move largest in column to beginning
745 //and do row part
746 {
747 CoinBigIndex i, k;
748
749 i = 0;
750 int iRow;
751 for ( iRow = 0; iRow < numberRows; iRow++ ) {
752 startRow[iRow] = i;
753 i += numberInRow[iRow];
754 }
755 CoinZeroN ( numberInRow, numberRows );
756 int iColumn;
757 for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
758 int number = numberInColumn[iColumn];
759
760 if ( number ) {
761 CoinBigIndex first = startColumn[iColumn];
762 CoinBigIndex largest = first;
763 int iRowSave = indexRow[first];
764 CoinFactorizationDouble valueSave = element[first];
765 double valueLargest = fabs ( valueSave );
766 int iLook = numberInRow[iRowSave];
767
768 numberInRow[iRowSave] = iLook + 1;
769 indexColumn[startRow[iRowSave] + iLook] = iColumn;
770 for ( k = first + 1; k < first + number; k++ ) {
771 int iRow = indexRow[k];
772 int iLook = numberInRow[iRow];
773
774 numberInRow[iRow] = iLook + 1;
775 indexColumn[startRow[iRow] + iLook] = iColumn;
776 CoinFactorizationDouble value = element[k];
777 double valueAbs = fabs ( value );
778
779 if ( valueAbs > valueLargest ) {
780 valueLargest = valueAbs;
781 largest = k;
782 }
783 }
784 indexRow[first] = indexRow[largest];
785 element[first] = element[largest];
786 indexRow[largest] = iRowSave;
787 element[largest] = valueSave;
788 }
789 }
790 }
791 case 3: //links and initialize pivots
792 {
793 int *lastRow = lastRow_.array();
794 int *nextRow = nextRow_.array();
795 int *lastColumn = lastColumn_.array();
796 int *nextColumn = nextColumn_.array();
797
798 CoinFillN ( firstCount_.array(), biggerDimension_ + 2, -1 );
799 CoinFillN ( pivotColumn_.array(), numberColumns_, -1 );
800 CoinZeroN ( numberInColumnPlus, maximumColumnsExtra_ + 1 );
801 int iRow;
802 for ( iRow = 0; iRow < numberRows; iRow++ ) {
803 lastRow[iRow] = iRow - 1;
804 nextRow[iRow] = iRow + 1;
805 int number = numberInRow[iRow];
806
807 addLink ( iRow, number );
808 }
809 lastRow[maximumRowsExtra_] = numberRows - 1;
810 nextRow[maximumRowsExtra_] = 0;
811 lastRow[0] = maximumRowsExtra_;
812 nextRow[numberRows - 1] = maximumRowsExtra_;
813 startRow[maximumRowsExtra_] = numberElements;
814 int iColumn;
815 for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
816 lastColumn[iColumn] = iColumn - 1;
817 nextColumn[iColumn] = iColumn + 1;
818 int number = numberInColumn[iColumn];
819
820 addLink ( iColumn + numberRows, number );
821 }
822 lastColumn[maximumColumnsExtra_] = numberColumns - 1;
823 nextColumn[maximumColumnsExtra_] = 0;
824 lastColumn[0] = maximumColumnsExtra_;
825 if (numberColumns)
826 nextColumn[numberColumns - 1] = maximumColumnsExtra_;
827 startColumn[maximumColumnsExtra_] = numberElements;
828 }
829 break;
830 case 4: //move largest in column to beginning
831 {
832 CoinBigIndex i, k;
833 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
834 int iColumn;
835 int iRow;
836 for ( iRow = 0; iRow < numberRows; iRow++ ) {
837 if( numberInRow[iRow]>=0) {
838 // zero count
839 numberInRow[iRow]=0;
840 } else {
841 // empty
842 //numberInRow[iRow]=-1; already that
843 }
844 }
845 //CoinZeroN ( numberInColumnPlus, maximumColumnsExtra_ + 1 );
846 for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
847 int number = numberInColumn[iColumn];
848
849 if ( number ) {
850 // use pivotRegion and startRow for remaining elements
851 CoinBigIndex first = startColumn[iColumn];
852 CoinBigIndex largest = -1;
853
854 double valueLargest = -1.0;
855 int nOther=0;
856 k = first;
857 CoinBigIndex end = first+number;
858 for ( ; k < end; k++ ) {
859 int iRow = indexRow[k];
860 assert (iRow<numberRows_);
861 CoinFactorizationDouble value = element[k];
862 if (numberInRow[iRow]>=0) {
863 numberInRow[iRow]++;
864 double valueAbs = fabs ( value );
865 if ( valueAbs > valueLargest ) {
866 valueLargest = valueAbs;
867 largest = nOther;
868 }
869 startRow[nOther]=iRow;
870 pivotRegion[nOther++]=value;
871 } else {
872 indexRow[first] = iRow;
873 element[first++] = value;
874 }
875 }
876 numberInColumnPlus[iColumn]=first-startColumn[iColumn];
877 startColumn[iColumn]=first;
878 //largest
879 if (largest>=0) {
880 indexRow[first] = startRow[largest];
881 element[first++] = pivotRegion[largest];
882 }
883 for (k=0;k<nOther;k++) {
884 if (k!=largest) {
885 indexRow[first] = startRow[k];
886 element[first++] = pivotRegion[k];
887 }
888 }
889 numberInColumn[iColumn]=first-startColumn[iColumn];
890 }
891 }
892 //and do row part
893 i = 0;
894 for ( iRow = 0; iRow < numberRows; iRow++ ) {
895 startRow[iRow] = i;
896 int n=numberInRow[iRow];
897 if (n>0) {
898 numberInRow[iRow]=0;
899 i += n;
900 }
901 }
902 for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
903 int number = numberInColumn[iColumn];
904
905 if ( number ) {
906 CoinBigIndex first = startColumn[iColumn];
907 for ( k = first; k < first + number; k++ ) {
908 int iRow = indexRow[k];
909 int iLook = numberInRow[iRow];
910
911 numberInRow[iRow] = iLook + 1;
912 indexColumn[startRow[iRow] + iLook] = iColumn;
913 }
914 }
915 }
916 }
917 // modified 3
918 {
919 //set markRow so no rows updated
920 //CoinFillN ( markRow_.array(), numberRows_, -1 );
921 int *lastColumn = lastColumn_.array();
922 int *nextColumn = nextColumn_.array();
923 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
924
925 int iRow;
926 int numberGood=0;
927 startColumnL_.array()[0] = 0; //for luck
928 for ( iRow = 0; iRow < numberRows; iRow++ ) {
929 int number = numberInRow[iRow];
930 if (number<0) {
931 numberInRow[iRow]=0;
932 pivotRegion[numberGood++]=slackValue_;
933 }
934 }
935 int iColumn;
936 for ( iColumn = 0 ; iColumn < numberColumns_; iColumn++ ) {
937 lastColumn[iColumn] = iColumn - 1;
938 nextColumn[iColumn] = iColumn + 1;
939 int number = numberInColumn[iColumn];
940 deleteLink(iColumn+numberRows);
941 addLink ( iColumn + numberRows, number );
942 }
943 lastColumn[maximumColumnsExtra_] = numberColumns - 1;
944 nextColumn[maximumColumnsExtra_] = 0;
945 lastColumn[0] = maximumColumnsExtra_;
946 if (numberColumns)
947 nextColumn[numberColumns - 1] = maximumColumnsExtra_;
948 startColumn[maximumColumnsExtra_] = numberElements;
949 }
950 } /* endswitch */
951}
952
953//Does most of factorization
954int
955CoinFactorization::factor ( )
956{
957 int * lastColumn = lastColumn_.array();
958 int * lastRow = lastRow_.array();
959 //sparse
960 status_ = factorSparse ( );
961 switch ( status_ ) {
962 case 0: //finished
963 totalElements_ = 0;
964 {
965 int * pivotColumn = pivotColumn_.array();
966 if ( numberGoodU_ < numberRows_ ) {
967 int i, k;
968 // Clean out unset nextRow
969 int * nextRow = nextRow_.array();
970 //int nSing =0;
971 k=nextRow[maximumRowsExtra_];
972 while (k!=maximumRowsExtra_) {
973 int iRow = k;
974 k=nextRow[k];
975 //nSing++;
976 nextRow[iRow]=-1;
977 }
978 //assert (nSing);
979 //printf("%d singularities - good %d rows %d\n",nSing,
980 // numberGoodU_,numberRows_);
981 // Now nextRow has -1 or sequence into numberGoodU_;
982 int * permuteA = permute_.array();
983 for ( i = 0; i < numberRows_; i++ ) {
984 int iGood= nextRow[i];
985 if (iGood>=0)
986 permuteA[iGood]=i;
987 }
988
989 // swap arrays
990 permute_.swap(nextRow_);
991 int * permute = permute_.array();
992 for ( i = 0; i < numberRows_; i++ ) {
993 lastRow[i] = -1;
994 }
995 for ( i = 0; i < numberColumns_; i++ ) {
996 lastColumn[i] = -1;
997 }
998 for ( i = 0; i < numberGoodU_; i++ ) {
999 int goodRow = permuteA[i]; //valid pivot row
1000 int goodColumn = pivotColumn[i];
1001
1002 lastRow[goodRow] = goodColumn; //will now have -1 or column sequence
1003 lastColumn[goodColumn] = goodRow; //will now have -1 or row sequence
1004 }
1005 nextRow_.conditionalDelete();
1006 k = 0;
1007 //copy back and count
1008 for ( i = 0; i < numberRows_; i++ ) {
1009 permute[i] = lastRow[i];
1010 if ( permute[i] < 0 ) {
1011 //std::cout << i << " " <<permute[i] << std::endl;
1012 } else {
1013 k++;
1014 }
1015 }
1016 for ( i = 0; i < numberColumns_; i++ ) {
1017 pivotColumn[i] = lastColumn[i];
1018 }
1019 if ((messageLevel_&4)!=0)
1020 std::cout <<"Factorization has "<<numberRows_-k
1021 <<" singularities"<<std::endl;
1022 status_ = -1;
1023 }
1024 }
1025 break;
1026 // dense
1027 case 2:
1028 status_=factorDense();
1029 if(!status_)
1030 break;
1031 default:
1032 //singular ? or some error
1033 if ((messageLevel_&4)!=0)
1034 std::cout << "Error " << status_ << std::endl;
1035 break;
1036 } /* endswitch */
1037 //clean up
1038 if ( !status_ ) {
1039 if ( (messageLevel_ & 16)&&numberCompressions_)
1040 std::cout<<" Factorization did "<<numberCompressions_
1041 <<" compressions"<<std::endl;
1042 if ( numberCompressions_ > 10 ) {
1043 areaFactor_ *= 1.1;
1044 }
1045 numberCompressions_=0;
1046 cleanup ( );
1047 }
1048 return status_;
1049}
1050
1051
1052// pivotRowSingleton. Does one pivot on Row Singleton in factorization
1053bool
1054CoinFactorization::pivotRowSingleton ( int pivotRow,
1055 int pivotColumn )
1056{
1057 //store pivot columns (so can easily compress)
1058 CoinBigIndex * startColumnU = startColumnU_.array();
1059 CoinBigIndex startColumn = startColumnU[pivotColumn];
1060 int * numberInRow = numberInRow_.array();
1061 int * numberInColumn = numberInColumn_.array();
1062 int numberDoColumn = numberInColumn[pivotColumn] - 1;
1063 CoinBigIndex endColumn = startColumn + numberDoColumn + 1;
1064 CoinBigIndex pivotRowPosition = startColumn;
1065 int * indexRowU = indexRowU_.array();
1066 int iRow = indexRowU[pivotRowPosition];
1067 CoinBigIndex * startRowU = startRowU_.array();
1068 int * nextRow = nextRow_.array();
1069 int * lastRow = lastRow_.array();
1070
1071 while ( iRow != pivotRow ) {
1072 pivotRowPosition++;
1073 iRow = indexRowU[pivotRowPosition];
1074 } /* endwhile */
1075 assert ( pivotRowPosition < endColumn );
1076 //store column in L, compress in U and take column out
1077 CoinBigIndex l = lengthL_;
1078
1079 if ( l + numberDoColumn > lengthAreaL_ ) {
1080 //need more memory
1081 if ((messageLevel_&4)!=0)
1082 std::cout << "more memory needed in middle of invert" << std::endl;
1083 return false;
1084 }
1085 CoinBigIndex * startColumnL = startColumnL_.array();
1086 CoinFactorizationDouble * elementL = elementL_.array();
1087 int * indexRowL = indexRowL_.array();
1088 startColumnL[numberGoodL_] = l; //for luck and first time
1089 numberGoodL_++;
1090 startColumnL[numberGoodL_] = l + numberDoColumn;
1091 lengthL_ += numberDoColumn;
1092 CoinFactorizationDouble *elementU = elementU_.array();
1093 CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
1094 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1095
1096 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1097 CoinBigIndex i;
1098
1099 int * indexColumnU = indexColumnU_.array();
1100 for ( i = startColumn; i < pivotRowPosition; i++ ) {
1101 int iRow = indexRowU[i];
1102
1103 indexRowL[l] = iRow;
1104 elementL[l] = elementU[i] * pivotMultiplier;
1105 l++;
1106 //take out of row list
1107 CoinBigIndex start = startRowU[iRow];
1108 int iNumberInRow = numberInRow[iRow];
1109 CoinBigIndex end = start + iNumberInRow;
1110 CoinBigIndex where = start;
1111
1112 while ( indexColumnU[where] != pivotColumn ) {
1113 where++;
1114 } /* endwhile */
1115 assert ( where < end );
1116 indexColumnU[where] = indexColumnU[end - 1];
1117 iNumberInRow--;
1118 numberInRow[iRow] = iNumberInRow;
1119 deleteLink ( iRow );
1120 addLink ( iRow, iNumberInRow );
1121 }
1122 for ( i = pivotRowPosition + 1; i < endColumn; i++ ) {
1123 int iRow = indexRowU[i];
1124
1125 indexRowL[l] = iRow;
1126 elementL[l] = elementU[i] * pivotMultiplier;
1127 l++;
1128 //take out of row list
1129 CoinBigIndex start = startRowU[iRow];
1130 int iNumberInRow = numberInRow[iRow];
1131 CoinBigIndex end = start + iNumberInRow;
1132 CoinBigIndex where = start;
1133
1134 while ( indexColumnU[where] != pivotColumn ) {
1135 where++;
1136 } /* endwhile */
1137 assert ( where < end );
1138 indexColumnU[where] = indexColumnU[end - 1];
1139 iNumberInRow--;
1140 numberInRow[iRow] = iNumberInRow;
1141 deleteLink ( iRow );
1142 addLink ( iRow, iNumberInRow );
1143 }
1144 numberInColumn[pivotColumn] = 0;
1145 //modify linked list for pivots
1146 numberInRow[pivotRow] = 0;
1147 deleteLink ( pivotRow );
1148 deleteLink ( pivotColumn + numberRows_ );
1149 //take out this bit of indexColumnU
1150 int next = nextRow[pivotRow];
1151 int last = lastRow[pivotRow];
1152
1153 nextRow[last] = next;
1154 lastRow[next] = last;
1155 lastRow[pivotRow] =-2;
1156 nextRow[pivotRow] = numberGoodU_; //use for permute
1157 return true;
1158}
1159
1160// pivotColumnSingleton. Does one pivot on Column Singleton in factorization
1161bool
1162CoinFactorization::pivotColumnSingleton ( int pivotRow,
1163 int pivotColumn )
1164{
1165 int * numberInRow = numberInRow_.array();
1166 int * numberInColumn = numberInColumn_.array();
1167 int * numberInColumnPlus = numberInColumnPlus_.array();
1168 //store pivot columns (so can easily compress)
1169 int numberDoRow = numberInRow[pivotRow] - 1;
1170 CoinBigIndex * startColumnU = startColumnU_.array();
1171 CoinBigIndex startColumn = startColumnU[pivotColumn];
1172 int put = 0;
1173 CoinBigIndex * startRowU = startRowU_.array();
1174 CoinBigIndex startRow = startRowU[pivotRow];
1175 CoinBigIndex endRow = startRow + numberDoRow + 1;
1176 int * indexColumnU = indexColumnU_.array();
1177 int * indexRowU = indexRowU_.array();
1178 int * saveColumn = saveColumn_.array();
1179 CoinBigIndex i;
1180
1181 for ( i = startRow; i < endRow; i++ ) {
1182 int iColumn = indexColumnU[i];
1183
1184 if ( iColumn != pivotColumn ) {
1185 saveColumn[put++] = iColumn;
1186 }
1187 }
1188 int * nextRow = nextRow_.array();
1189 int * lastRow = lastRow_.array();
1190 //take out this bit of indexColumnU
1191 int next = nextRow[pivotRow];
1192 int last = lastRow[pivotRow];
1193
1194 nextRow[last] = next;
1195 lastRow[next] = last;
1196 nextRow[pivotRow] = numberGoodU_; //use for permute
1197 lastRow[pivotRow] =-2; //mark
1198 //clean up counts
1199 CoinFactorizationDouble *elementU = elementU_.array();
1200 CoinFactorizationDouble pivotElement = elementU[startColumn];
1201
1202 pivotRegion_.array()[numberGoodU_] = 1.0 / pivotElement;
1203 numberInColumn[pivotColumn] = 0;
1204 //totalElements_ --;
1205 //numberInColumnPlus[pivotColumn]++;
1206 //move pivot row in other columns to safe zone
1207 for ( i = 0; i < numberDoRow; i++ ) {
1208 int iColumn = saveColumn[i];
1209
1210 if ( numberInColumn[iColumn] ) {
1211 int number = numberInColumn[iColumn] - 1;
1212
1213 //modify linked list
1214 deleteLink ( iColumn + numberRows_ );
1215 addLink ( iColumn + numberRows_, number );
1216 //move pivot row element
1217 if ( number ) {
1218 CoinBigIndex start = startColumnU[iColumn];
1219 CoinBigIndex pivot = start;
1220 int iRow = indexRowU[pivot];
1221 while ( iRow != pivotRow ) {
1222 pivot++;
1223 iRow = indexRowU[pivot];
1224 }
1225 assert (pivot < startColumnU[iColumn] +
1226 numberInColumn[iColumn]);
1227 if ( pivot != start ) {
1228 //move largest one up
1229 CoinFactorizationDouble value = elementU[start];
1230
1231 iRow = indexRowU[start];
1232 elementU[start] = elementU[pivot];
1233 indexRowU[start] = indexRowU[pivot];
1234 elementU[pivot] = elementU[start + 1];
1235 indexRowU[pivot] = indexRowU[start + 1];
1236 elementU[start + 1] = value;
1237 indexRowU[start + 1] = iRow;
1238 } else {
1239 //find new largest element
1240 int iRowSave = indexRowU[start + 1];
1241 CoinFactorizationDouble valueSave = elementU[start + 1];
1242 double valueLargest = fabs ( valueSave );
1243 CoinBigIndex end = start + numberInColumn[iColumn];
1244 CoinBigIndex largest = start + 1;
1245
1246 CoinBigIndex k;
1247 for ( k = start + 2; k < end; k++ ) {
1248 CoinFactorizationDouble value = elementU[k];
1249 double valueAbs = fabs ( value );
1250
1251 if ( valueAbs > valueLargest ) {
1252 valueLargest = valueAbs;
1253 largest = k;
1254 }
1255 }
1256 indexRowU[start + 1] = indexRowU[largest];
1257 elementU[start + 1] = elementU[largest];
1258 indexRowU[largest] = iRowSave;
1259 elementU[largest] = valueSave;
1260 }
1261 }
1262 //clean up counts
1263 numberInColumn[iColumn]--;
1264 numberInColumnPlus[iColumn]++;
1265 startColumnU[iColumn]++;
1266 //totalElements_--;
1267 }
1268 }
1269 //modify linked list for pivots
1270 deleteLink ( pivotRow );
1271 deleteLink ( pivotColumn + numberRows_ );
1272 numberInRow[pivotRow] = 0;
1273 //put in dummy pivot in L
1274 CoinBigIndex l = lengthL_;
1275
1276 CoinBigIndex * startColumnL = startColumnL_.array();
1277 startColumnL[numberGoodL_] = l; //for luck and first time
1278 numberGoodL_++;
1279 startColumnL[numberGoodL_] = l;
1280 return true;
1281}
1282
1283
1284// getColumnSpace. Gets space for one Column with given length
1285//may have to do compression (returns true)
1286//also moves existing vector
1287bool
1288CoinFactorization::getColumnSpace ( int iColumn,
1289 int extraNeeded )
1290{
1291 int * numberInColumn = numberInColumn_.array();
1292 int * numberInColumnPlus = numberInColumnPlus_.array();
1293 int * nextColumn = nextColumn_.array();
1294 int * lastColumn = lastColumn_.array();
1295 int number = numberInColumnPlus[iColumn] +
1296 numberInColumn[iColumn];
1297 CoinBigIndex * startColumnU = startColumnU_.array();
1298 CoinBigIndex space = lengthAreaU_ - startColumnU[maximumColumnsExtra_];
1299 CoinFactorizationDouble *elementU = elementU_.array();
1300 int * indexRowU = indexRowU_.array();
1301
1302 if ( space < extraNeeded + number + 4 ) {
1303 //compression
1304 int iColumn = nextColumn[maximumColumnsExtra_];
1305 CoinBigIndex put = 0;
1306
1307 while ( iColumn != maximumColumnsExtra_ ) {
1308 //move
1309 CoinBigIndex get;
1310 CoinBigIndex getEnd;
1311
1312 if ( startColumnU[iColumn] >= 0 ) {
1313 get = startColumnU[iColumn]
1314 - numberInColumnPlus[iColumn];
1315 getEnd = startColumnU[iColumn] + numberInColumn[iColumn];
1316 startColumnU[iColumn] = put + numberInColumnPlus[iColumn];
1317 } else {
1318 get = -startColumnU[iColumn];
1319 getEnd = get + numberInColumn[iColumn];
1320 startColumnU[iColumn] = -put;
1321 }
1322 CoinBigIndex i;
1323 for ( i = get; i < getEnd; i++ ) {
1324 indexRowU[put] = indexRowU[i];
1325 elementU[put] = elementU[i];
1326 put++;
1327 }
1328 iColumn = nextColumn[iColumn];
1329 } /* endwhile */
1330 numberCompressions_++;
1331 startColumnU[maximumColumnsExtra_] = put;
1332 space = lengthAreaU_ - put;
1333 if ( extraNeeded == COIN_INT_MAX >> 1 ) {
1334 return true;
1335 }
1336 if ( space < extraNeeded + number + 2 ) {
1337 //need more space
1338 //if we can allocate bigger then do so and copy
1339 //if not then return so code can start again
1340 status_ = -99;
1341 return false;
1342 }
1343 }
1344 CoinBigIndex put = startColumnU[maximumColumnsExtra_];
1345 int next = nextColumn[iColumn];
1346 int last = lastColumn[iColumn];
1347
1348 if ( extraNeeded || next != maximumColumnsExtra_ ) {
1349 //out
1350 nextColumn[last] = next;
1351 lastColumn[next] = last;
1352 //in at end
1353 last = lastColumn[maximumColumnsExtra_];
1354 nextColumn[last] = iColumn;
1355 lastColumn[maximumColumnsExtra_] = iColumn;
1356 lastColumn[iColumn] = last;
1357 nextColumn[iColumn] = maximumColumnsExtra_;
1358 //move
1359 CoinBigIndex get = startColumnU[iColumn]
1360 - numberInColumnPlus[iColumn];
1361
1362 startColumnU[iColumn] = put + numberInColumnPlus[iColumn];
1363 if ( number < 50 ) {
1364 int *indexRow = indexRowU;
1365 CoinFactorizationDouble *element = elementU;
1366 int i = 0;
1367
1368 if ( ( number & 1 ) != 0 ) {
1369 element[put] = element[get];
1370 indexRow[put] = indexRow[get];
1371 i = 1;
1372 }
1373 for ( ; i < number; i += 2 ) {
1374 CoinFactorizationDouble value0 = element[get + i];
1375 CoinFactorizationDouble value1 = element[get + i + 1];
1376 int index0 = indexRow[get + i];
1377 int index1 = indexRow[get + i + 1];
1378
1379 element[put + i] = value0;
1380 element[put + i + 1] = value1;
1381 indexRow[put + i] = index0;
1382 indexRow[put + i + 1] = index1;
1383 }
1384 } else {
1385 CoinMemcpyN ( &indexRowU[get], number, &indexRowU[put] );
1386 CoinMemcpyN ( &elementU[get], number, &elementU[put] );
1387 }
1388 put += number;
1389 get += number;
1390 //add 2 for luck
1391 startColumnU[maximumColumnsExtra_] = put + extraNeeded + 2;
1392 if (startColumnU[maximumColumnsExtra_]>lengthAreaU_) {
1393 // get more memory
1394#ifdef CLP_DEVELOP
1395 printf("put %d, needed %d, start %d, length %d\n",
1396 put,extraNeeded,startColumnU[maximumColumnsExtra_],
1397 lengthAreaU_);
1398#endif
1399 return false;
1400 }
1401 } else {
1402 //take off space
1403 startColumnU[maximumColumnsExtra_] = startColumnU[last] +
1404 numberInColumn[last];
1405 }
1406 return true;
1407}
1408
1409// getRowSpace. Gets space for one Row with given length
1410//may have to do compression (returns true)
1411//also moves existing vector
1412bool
1413CoinFactorization::getRowSpace ( int iRow,
1414 int extraNeeded )
1415{
1416 int * numberInRow = numberInRow_.array();
1417 int number = numberInRow[iRow];
1418 CoinBigIndex * startRowU = startRowU_.array();
1419 CoinBigIndex space = lengthAreaU_ - startRowU[maximumRowsExtra_];
1420 int * nextRow = nextRow_.array();
1421 int * lastRow = lastRow_.array();
1422 int * indexColumnU = indexColumnU_.array();
1423
1424 if ( space < extraNeeded + number + 2 ) {
1425 //compression
1426 int iRow = nextRow[maximumRowsExtra_];
1427 CoinBigIndex put = 0;
1428
1429 while ( iRow != maximumRowsExtra_ ) {
1430 //move
1431 CoinBigIndex get = startRowU[iRow];
1432 CoinBigIndex getEnd = startRowU[iRow] + numberInRow[iRow];
1433
1434 startRowU[iRow] = put;
1435 CoinBigIndex i;
1436 for ( i = get; i < getEnd; i++ ) {
1437 indexColumnU[put] = indexColumnU[i];
1438 put++;
1439 }
1440 iRow = nextRow[iRow];
1441 } /* endwhile */
1442 numberCompressions_++;
1443 startRowU[maximumRowsExtra_] = put;
1444 space = lengthAreaU_ - put;
1445 if ( space < extraNeeded + number + 2 ) {
1446 //need more space
1447 //if we can allocate bigger then do so and copy
1448 //if not then return so code can start again
1449 status_ = -99;
1450 return false;
1451 }
1452 }
1453 CoinBigIndex put = startRowU[maximumRowsExtra_];
1454 int next = nextRow[iRow];
1455 int last = lastRow[iRow];
1456
1457 //out
1458 nextRow[last] = next;
1459 lastRow[next] = last;
1460 //in at end
1461 last = lastRow[maximumRowsExtra_];
1462 nextRow[last] = iRow;
1463 lastRow[maximumRowsExtra_] = iRow;
1464 lastRow[iRow] = last;
1465 nextRow[iRow] = maximumRowsExtra_;
1466 //move
1467 CoinBigIndex get = startRowU[iRow];
1468
1469 startRowU[iRow] = put;
1470 while ( number ) {
1471 number--;
1472 indexColumnU[put] = indexColumnU[get];
1473 put++;
1474 get++;
1475 } /* endwhile */
1476 //add 4 for luck
1477 startRowU[maximumRowsExtra_] = put + extraNeeded + 4;
1478 return true;
1479}
1480
1481#if COIN_ONE_ETA_COPY
1482/* Reorders U so contiguous and in order (if there is space)
1483 Returns true if it could */
1484bool
1485CoinFactorization::reorderU()
1486{
1487#if 1
1488 return false;
1489#else
1490 if (numberRows_!=numberColumns_)
1491 return false;
1492 CoinBigIndex * startColumnU = startColumnU_.array();
1493 int * numberInColumn = numberInColumn_.array();
1494 int * numberInColumnPlus = numberInColumnPlus_.array();
1495 int iColumn;
1496 CoinBigIndex put = 0;
1497 for (iColumn =0;iColumn<numberRows_;iColumn++)
1498 put += numberInColumnPlus[iColumn];
1499 CoinBigIndex space = lengthAreaU_ - startColumnU[maximumColumnsExtra_];
1500 if (space<put) {
1501 //printf("Space %d out of %d - needed %d\n",
1502 // space,lengthAreaU_,put);
1503 return false;
1504 }
1505 int *indexRowU = indexRowU_.array();
1506 CoinFactorizationDouble *elementU = elementU_.array();
1507 int * pivotColumn = pivotColumn_.array();
1508 put = startColumnU[maximumColumnsExtra_];
1509 for (int jColumn =0;jColumn<numberRows_;jColumn++) {
1510 iColumn = pivotColumn[jColumn];
1511 int n = numberInColumnPlus[iColumn];
1512 CoinBigIndex getEnd = startColumnU[iColumn];
1513 CoinBigIndex get = getEnd - n;
1514 startColumnU[iColumn] = put;
1515 numberInColumn[jColumn]=n;
1516 CoinBigIndex i;
1517 for ( i = get; i < getEnd; i++ ) {
1518 indexRowU[put] = indexRowU[i];
1519 elementU[put] = elementU[i];
1520 put++;
1521 }
1522 }
1523 // and pack down
1524 put = 0;
1525 for (int jColumn =0;jColumn<numberRows_;jColumn++) {
1526 iColumn = pivotColumn[jColumn];
1527 int n = numberInColumn[jColumn];
1528 CoinBigIndex get = startColumnU[iColumn];
1529 CoinBigIndex getEnd = get+n;
1530 CoinBigIndex i;
1531 for ( i = get; i < getEnd; i++ ) {
1532 indexRowU[put] = indexRowU[i];
1533 elementU[put] = elementU[i];
1534 put++;
1535 }
1536 }
1537 put=0;
1538 for (iColumn =0;iColumn<numberRows_;iColumn++) {
1539 int n = numberInColumn[iColumn];
1540 startColumnU[iColumn]=put;
1541 put += n;
1542 //numberInColumnPlus[iColumn]=n;
1543 //numberInColumn[iColumn]=0; // necessary?
1544 //pivotColumn[iColumn]=iColumn;
1545 }
1546#if 0
1547 // reset nextColumn - probably not necessary
1548 int * nextColumn = nextColumn_.array();
1549 nextColumn[maximumColumnsExtra_]=0;
1550 nextColumn[numberRows_-1] = maximumColumnsExtra_;
1551 for (iColumn=0;iColumn<numberRows_-1;iColumn++)
1552 nextColumn[iColumn]=iColumn+1;
1553 // swap arrays
1554 numberInColumn_.swap(numberInColumnPlus_);
1555#endif
1556 //return false;
1557 return true;
1558#endif
1559}
1560#endif
1561// cleanup. End of factorization
1562void
1563CoinFactorization::cleanup ( )
1564{
1565#if COIN_ONE_ETA_COPY
1566 bool compressDone = reorderU();
1567 if (!compressDone) {
1568 getColumnSpace ( 0, COIN_INT_MAX >> 1 ); //compress
1569 // swap arrays
1570 numberInColumn_.swap(numberInColumnPlus_);
1571 }
1572#else
1573 getColumnSpace ( 0, COIN_INT_MAX >> 1 ); //compress
1574 // swap arrays
1575 numberInColumn_.swap(numberInColumnPlus_);
1576#endif
1577 CoinBigIndex * startColumnU = startColumnU_.array();
1578 CoinBigIndex lastU = startColumnU[maximumColumnsExtra_];
1579
1580 //free some memory here
1581 saveColumn_.conditionalDelete();
1582 markRow_.conditionalDelete() ;
1583 //firstCount_.conditionalDelete() ;
1584 nextCount_.conditionalDelete() ;
1585 lastCount_.conditionalDelete() ;
1586 int * numberInRow = numberInRow_.array();
1587 int * numberInColumn = numberInColumn_.array();
1588 int * numberInColumnPlus = numberInColumnPlus_.array();
1589 //make column starts OK
1590 //for best cache behavior get in order (last pivot at bottom of space)
1591 //that will need thinking about
1592 //use nextRow for permutation (as that is what it is)
1593 int i;
1594
1595#ifndef NDEBUG
1596 {
1597 if (numberGoodU_<numberRows_)
1598 abort();
1599 char * mark = new char[numberRows_];
1600 memset(mark,0,numberRows_);
1601 int * array;
1602 array = nextRow_.array();
1603 for ( i = 0; i < numberRows_; i++ ) {
1604 int k = array[i];
1605 if(k<0||k>=numberRows_)
1606 printf("Bad a %d %d\n",i,k);
1607 assert(k>=0&&k<numberRows_);
1608 if(mark[k]==1)
1609 printf("Bad a %d %d\n",i,k);
1610 mark[k]=1;
1611 }
1612 for ( i = 0; i < numberRows_; i++ ) {
1613 assert(mark[i]==1);
1614 if(mark[i]!=1)
1615 printf("Bad b %d\n",i);
1616 }
1617 delete [] mark;
1618 }
1619#endif
1620 // swap arrays
1621 permute_.swap(nextRow_);
1622 //safety feature
1623 int * permute = permute_.array();
1624 permute[numberRows_] = 0;
1625 permuteBack_.conditionalNew( maximumRowsExtra_ + 1);
1626 int * permuteBack = permuteBack_.array();
1627#ifdef ZEROFAULT
1628 memset(permuteBack_.array(),'w',(maximumRowsExtra_+1)*sizeof(int));
1629#endif
1630 for ( i = 0; i < numberRows_; i++ ) {
1631 int iRow = permute[i];
1632
1633 permuteBack[iRow] = i;
1634 }
1635 //redo nextRow_
1636
1637#ifndef NDEBUG
1638 for ( i = 0; i < numberRows_; i++ ) {
1639 assert (permute[i]>=0&&permute[i]<numberRows_);
1640 assert (permuteBack[i]>=0&&permuteBack[i]<numberRows_);
1641 }
1642#endif
1643#if COIN_ONE_ETA_COPY
1644 if (!compressDone) {
1645#endif
1646 // Redo total elements
1647 totalElements_=0;
1648 for ( i = 0; i < numberColumns_; i++ ) {
1649 int number = numberInColumn[i];
1650 totalElements_ += number;
1651 startColumnU[i] -= number;
1652 }
1653#if COIN_ONE_ETA_COPY
1654 }
1655#endif
1656 int numberU = 0;
1657
1658 pivotColumnBack_.conditionalNew( maximumRowsExtra_ + 1);
1659#ifdef ZEROFAULT
1660 memset(pivotColumnBack(),'q',(maximumRowsExtra_+1)*sizeof(int));
1661#endif
1662 const int * pivotColumn = pivotColumn_.array();
1663 int * pivotColumnB = pivotColumnBack();
1664 int *indexColumnU = indexColumnU_.array();
1665 CoinBigIndex *startColumn = startColumnU;
1666 int *indexRowU = indexRowU_.array();
1667 CoinFactorizationDouble *elementU = elementU_.array();
1668#if COIN_ONE_ETA_COPY
1669 if (!compressDone) {
1670#endif
1671 for ( i = 0; i < numberColumns_; i++ ) {
1672 int iColumn = pivotColumn[i];
1673
1674 pivotColumnB[iColumn] = i;
1675 if ( iColumn >= 0 ) {
1676 //wanted
1677 if ( numberU != iColumn ) {
1678 numberInColumnPlus[iColumn] = numberU;
1679 } else {
1680 numberInColumnPlus[iColumn] = -1; //already in correct place
1681 }
1682 numberU++;
1683 }
1684 }
1685 for ( i = 0; i < numberColumns_; i++ ) {
1686 int number = numberInColumn[i]; //always 0?
1687 int where = numberInColumnPlus[i];
1688
1689 numberInColumnPlus[i] = -1;
1690 CoinBigIndex start = startColumnU[i];
1691
1692 while ( where >= 0 ) {
1693 //put where it should be
1694 int numberNext = numberInColumn[where]; //always 0?
1695 int whereNext = numberInColumnPlus[where];
1696 CoinBigIndex startNext = startColumnU[where];
1697
1698 numberInColumn[where] = number;
1699 numberInColumnPlus[where] = -1;
1700 startColumnU[where] = start;
1701 number = numberNext;
1702 where = whereNext;
1703 start = startNext;
1704 } /* endwhile */
1705 }
1706 //sort - using indexColumn
1707 CoinFillN ( indexColumnU_.array(), lastU, -1 );
1708 CoinBigIndex k = 0;
1709
1710 for ( i = numberSlacks_; i < numberRows_; i++ ) {
1711 CoinBigIndex start = startColumn[i];
1712 CoinBigIndex end = start + numberInColumn[i];
1713
1714 CoinBigIndex j;
1715 for ( j = start; j < end; j++ ) {
1716 indexColumnU[j] = k++;
1717 }
1718 }
1719 for ( i = numberSlacks_; i < numberRows_; i++ ) {
1720 CoinBigIndex start = startColumn[i];
1721 CoinBigIndex end = start + numberInColumn[i];
1722
1723 CoinBigIndex j;
1724 for ( j = start; j < end; j++ ) {
1725 CoinBigIndex k = indexColumnU[j];
1726 int iRow = indexRowU[j];
1727 CoinFactorizationDouble element = elementU[j];
1728
1729 while ( k != -1 ) {
1730 CoinBigIndex kNext = indexColumnU[k];
1731 int iRowNext = indexRowU[k];
1732 CoinFactorizationDouble elementNext = elementU[k];
1733
1734 indexColumnU[k] = -1;
1735 indexRowU[k] = iRow;
1736 elementU[k] = element;
1737 k = kNext;
1738 iRow = iRowNext;
1739 element = elementNext;
1740 } /* endwhile */
1741 }
1742 }
1743 CoinZeroN ( startColumnU, numberSlacks_ );
1744 k = 0;
1745 for ( i = numberSlacks_; i < numberRows_; i++ ) {
1746 startColumnU[i] = k;
1747 k += numberInColumn[i];
1748 }
1749 maximumU_=k;
1750#if COIN_ONE_ETA_COPY
1751 } else {
1752 // U already OK
1753 for ( i = 0; i < numberColumns_; i++ ) {
1754 int iColumn = pivotColumn[i];
1755 pivotColumnB[iColumn] = i;
1756 }
1757 maximumU_=startColumnU[numberRows_-1]+
1758 numberInColumn[numberRows_-1];
1759 numberU=numberRows_;
1760 }
1761#endif
1762 if ( (messageLevel_ & 8)) {
1763 std::cout<<" length of U "<<totalElements_<<", length of L "<<lengthL_;
1764 if (numberDense_)
1765 std::cout<<" plus "<<numberDense_*numberDense_<<" from "<<numberDense_<<" dense rows";
1766 std::cout<<std::endl;
1767 }
1768 // and add L and dense
1769 totalElements_ += numberDense_*numberDense_+lengthL_;
1770 int * nextColumn = nextColumn_.array();
1771 int * lastColumn = lastColumn_.array();
1772 // See whether to have extra copy of R
1773 if (maximumU_>10*numberRows_||numberRows_<200) {
1774 // NO
1775 numberInColumnPlus_.conditionalDelete() ;
1776 } else {
1777 for ( i = 0; i < numberColumns_; i++ ) {
1778 lastColumn[i] = i - 1;
1779 nextColumn[i] = i + 1;
1780 numberInColumnPlus[i]=0;
1781 }
1782 nextColumn[numberColumns_ - 1] = maximumColumnsExtra_;
1783 lastColumn[maximumColumnsExtra_] = numberColumns_ - 1;
1784 nextColumn[maximumColumnsExtra_] = 0;
1785 lastColumn[0] = maximumColumnsExtra_;
1786 }
1787 numberU_ = numberU;
1788 numberGoodU_ = numberU;
1789 numberL_ = numberGoodL_;
1790#if COIN_DEBUG
1791 for ( i = 0; i < numberRows_; i++ ) {
1792 if ( permute[i] < 0 ) {
1793 std::cout << i << std::endl;
1794 abort ( );
1795 }
1796 }
1797#endif
1798 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
1799 for ( i = numberSlacks_; i < numberU; i++ ) {
1800 CoinBigIndex start = startColumnU[i];
1801 CoinBigIndex end = start + numberInColumn[i];
1802
1803 totalElements_ += numberInColumn[i];
1804 CoinBigIndex j;
1805
1806 for ( j = start; j < end; j++ ) {
1807 int iRow = indexRowU[j];
1808 iRow = permute[iRow];
1809 indexRowU[j] = iRow;
1810 numberInRow[iRow]++;
1811 }
1812 }
1813#if COIN_ONE_ETA_COPY
1814 if (numberRows_>=COIN_ONE_ETA_COPY) {
1815#endif
1816 //space for cross reference
1817 convertRowToColumnU_.conditionalNew( lengthAreaU_ );
1818 CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
1819 CoinBigIndex j = 0;
1820 CoinBigIndex *startRow = startRowU_.array();
1821
1822 int iRow;
1823 for ( iRow = 0; iRow < numberRows_; iRow++ ) {
1824 startRow[iRow] = j;
1825 j += numberInRow[iRow];
1826 }
1827 CoinBigIndex numberInU = j;
1828
1829 CoinZeroN ( numberInRow_.array(), numberRows_ );
1830
1831 for ( i = numberSlacks_; i < numberRows_; i++ ) {
1832 CoinBigIndex start = startColumnU[i];
1833 CoinBigIndex end = start + numberInColumn[i];
1834
1835 CoinFactorizationDouble pivotValue = pivotRegion[i];
1836
1837 CoinBigIndex j;
1838 for ( j = start; j < end; j++ ) {
1839 int iRow = indexRowU[j];
1840 int iLook = numberInRow[iRow];
1841
1842 numberInRow[iRow] = iLook + 1;
1843 CoinBigIndex k = startRow[iRow] + iLook;
1844
1845 indexColumnU[k] = i;
1846 convertRowToColumn[k] = j;
1847 //multiply by pivot
1848 elementU[j] *= pivotValue;
1849 }
1850 }
1851 int * nextRow = nextRow_.array();
1852 int * lastRow = lastRow_.array();
1853 for ( j = 0; j < numberRows_; j++ ) {
1854 lastRow[j] = j - 1;
1855 nextRow[j] = j + 1;
1856 }
1857 nextRow[numberRows_ - 1] = maximumRowsExtra_;
1858 lastRow[maximumRowsExtra_] = numberRows_ - 1;
1859 nextRow[maximumRowsExtra_] = 0;
1860 lastRow[0] = maximumRowsExtra_;
1861 startRow[maximumRowsExtra_] = numberInU;
1862#if COIN_ONE_ETA_COPY
1863 } else {
1864 // no row copy
1865 for ( i = numberSlacks_; i < numberU; i++ ) {
1866 CoinBigIndex start = startColumnU[i];
1867 CoinBigIndex end = start + numberInColumn[i];
1868
1869 CoinBigIndex j;
1870 CoinFactorizationDouble pivotValue = pivotRegion[i];
1871
1872 for ( j = start; j < end; j++ ) {
1873 //multiply by pivot
1874 elementU[j] *= pivotValue;
1875 }
1876 }
1877 }
1878#endif
1879
1880 int firstReal = numberRows_;
1881
1882 CoinBigIndex * startColumnL = startColumnL_.array();
1883 int * indexRowL = indexRowL_.array();
1884 for ( i = numberRows_ - 1; i >= 0; i-- ) {
1885 CoinBigIndex start = startColumnL[i];
1886 CoinBigIndex end = startColumnL[i + 1];
1887
1888 totalElements_ += end - start;
1889 if ( end > start ) {
1890 firstReal = i;
1891 CoinBigIndex j;
1892 for ( j = start; j < end; j++ ) {
1893 int iRow = indexRowL[j];
1894 iRow = permute[iRow];
1895 assert (iRow>firstReal);
1896 indexRowL[j] = iRow;
1897 }
1898 }
1899 }
1900 baseL_ = firstReal;
1901 numberL_ -= firstReal;
1902 factorElements_ = totalElements_;
1903 //can delete pivotRowL_ as not used
1904 pivotRowL_.conditionalDelete() ;
1905 //use L for R if room
1906 CoinBigIndex space = lengthAreaL_ - lengthL_;
1907 CoinBigIndex spaceUsed = lengthL_ + lengthU_;
1908
1909 int needed = ( spaceUsed + numberRows_ - 1 ) / numberRows_;
1910
1911 needed = needed * 2 * maximumPivots_;
1912 if ( needed < 2 * numberRows_ ) {
1913 needed = 2 * numberRows_;
1914 }
1915 if (numberInColumnPlus_.array()) {
1916 // Need double the space for R
1917 space = space/2;
1918 startColumnR_.conditionalNew( maximumPivots_ + 1 + maximumColumnsExtra_ + 1 );
1919 CoinBigIndex * startR = startColumnR_.array() + maximumPivots_+1;
1920 CoinZeroN (startR,(maximumColumnsExtra_+1));
1921 } else {
1922 startColumnR_.conditionalNew(maximumPivots_ + 1 );
1923 }
1924#ifdef ZEROFAULT
1925 memset(startColumnR_.array(),'z',(maximumPivots_ + 1)*sizeof(CoinBigIndex));
1926#endif
1927 if ( space >= needed ) {
1928 lengthR_ = 0;
1929 lengthAreaR_ = space;
1930 elementR_ = elementL_.array() + lengthL_;
1931 indexRowR_ = indexRowL_.array() + lengthL_;
1932 } else {
1933 lengthR_ = 0;
1934 lengthAreaR_ = space;
1935 elementR_ = elementL_.array() + lengthL_;
1936 indexRowR_ = indexRowL_.array() + lengthL_;
1937 if ((messageLevel_&4))
1938 std::cout<<"Factorization may need some increasing area space"
1939 <<std::endl;
1940 if ( areaFactor_ ) {
1941 areaFactor_ *= 1.1;
1942 } else {
1943 areaFactor_ = 1.1;
1944 }
1945 }
1946 numberR_ = 0;
1947}
1948// Returns areaFactor but adjusted for dense
1949double
1950CoinFactorization::adjustedAreaFactor() const
1951{
1952 double factor = areaFactor_;
1953 if (numberDense_&&areaFactor_>1.0) {
1954 double dense = numberDense_;
1955 dense *= dense;
1956 double withoutDense = totalElements_ - dense +1.0;
1957 factor *= 1.0 +dense/withoutDense;
1958 }
1959 return factor;
1960}
1961
1962// checkConsistency. Checks that row and column copies look OK
1963void
1964CoinFactorization::checkConsistency ( )
1965{
1966 bool bad = false;
1967
1968 int iRow;
1969 CoinBigIndex * startRowU = startRowU_.array();
1970 int * numberInRow = numberInRow_.array();
1971 int * numberInColumn = numberInColumn_.array();
1972 int * indexColumnU = indexColumnU_.array();
1973 int * indexRowU = indexRowU_.array();
1974 CoinBigIndex * startColumnU = startColumnU_.array();
1975 for ( iRow = 0; iRow < numberRows_; iRow++ ) {
1976 if ( numberInRow[iRow] ) {
1977 CoinBigIndex startRow = startRowU[iRow];
1978 CoinBigIndex endRow = startRow + numberInRow[iRow];
1979
1980 CoinBigIndex j;
1981 for ( j = startRow; j < endRow; j++ ) {
1982 int iColumn = indexColumnU[j];
1983 CoinBigIndex startColumn = startColumnU[iColumn];
1984 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
1985 bool found = false;
1986
1987 CoinBigIndex k;
1988 for ( k = startColumn; k < endColumn; k++ ) {
1989 if ( indexRowU[k] == iRow ) {
1990 found = true;
1991 break;
1992 }
1993 }
1994 if ( !found ) {
1995 bad = true;
1996 std::cout << "row " << iRow << " column " << iColumn << " Rows" << std::endl;
1997 }
1998 }
1999 }
2000 }
2001 int iColumn;
2002 for ( iColumn = 0; iColumn < numberColumns_; iColumn++ ) {
2003 if ( numberInColumn[iColumn] ) {
2004 CoinBigIndex startColumn = startColumnU[iColumn];
2005 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
2006
2007 CoinBigIndex j;
2008 for ( j = startColumn; j < endColumn; j++ ) {
2009 int iRow = indexRowU[j];
2010 CoinBigIndex startRow = startRowU[iRow];
2011 CoinBigIndex endRow = startRow + numberInRow[iRow];
2012 bool found = false;
2013
2014 CoinBigIndex k;
2015 for ( k = startRow; k < endRow; k++ ) {
2016 if ( indexColumnU[k] == iColumn ) {
2017 found = true;
2018 break;
2019 }
2020 }
2021 if ( !found ) {
2022 bad = true;
2023 std::cout << "row " << iRow << " column " << iColumn << " Columns" <<
2024 std::endl;
2025 }
2026 }
2027 }
2028 }
2029 if ( bad ) {
2030 abort ( );
2031 }
2032}
2033 // pivotOneOtherRow. When just one other row so faster
2034bool
2035CoinFactorization::pivotOneOtherRow ( int pivotRow,
2036 int pivotColumn )
2037{
2038 int * numberInRow = numberInRow_.array();
2039 int * numberInColumn = numberInColumn_.array();
2040 int * numberInColumnPlus = numberInColumnPlus_.array();
2041 int numberInPivotRow = numberInRow[pivotRow] - 1;
2042 CoinBigIndex * startRowU = startRowU_.array();
2043 CoinBigIndex * startColumnU = startColumnU_.array();
2044 CoinBigIndex startColumn = startColumnU[pivotColumn];
2045 CoinBigIndex startRow = startRowU[pivotRow];
2046 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
2047
2048 //take out this bit of indexColumnU
2049 int * nextRow = nextRow_.array();
2050 int * lastRow = lastRow_.array();
2051 int next = nextRow[pivotRow];
2052 int last = lastRow[pivotRow];
2053
2054 nextRow[last] = next;
2055 lastRow[next] = last;
2056 nextRow[pivotRow] = numberGoodU_; //use for permute
2057 lastRow[pivotRow] = -2;
2058 numberInRow[pivotRow] = 0;
2059 //store column in L, compress in U and take column out
2060 CoinBigIndex l = lengthL_;
2061
2062 if ( l + 1 > lengthAreaL_ ) {
2063 //need more memory
2064 if ((messageLevel_&4)!=0)
2065 std::cout << "more memory needed in middle of invert" << std::endl;
2066 return false;
2067 }
2068 //l+=currentAreaL_->elementByColumn-elementL_;
2069 //CoinBigIndex lSave=l;
2070 CoinBigIndex * startColumnL = startColumnL_.array();
2071 CoinFactorizationDouble * elementL = elementL_.array();
2072 int * indexRowL = indexRowL_.array();
2073 startColumnL[numberGoodL_] = l; //for luck and first time
2074 numberGoodL_++;
2075 startColumnL[numberGoodL_] = l + 1;
2076 lengthL_++;
2077 CoinFactorizationDouble pivotElement;
2078 CoinFactorizationDouble otherMultiplier;
2079 int otherRow;
2080 int * saveColumn = saveColumn_.array();
2081 CoinFactorizationDouble *elementU = elementU_.array();
2082 int * indexRowU = indexRowU_.array();
2083
2084 if ( indexRowU[startColumn] == pivotRow ) {
2085 pivotElement = elementU[startColumn];
2086 otherMultiplier = elementU[startColumn + 1];
2087 otherRow = indexRowU[startColumn + 1];
2088 } else {
2089 pivotElement = elementU[startColumn + 1];
2090 otherMultiplier = elementU[startColumn];
2091 otherRow = indexRowU[startColumn];
2092 }
2093 int numberSave = numberInRow[otherRow];
2094 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
2095
2096 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
2097 pivotRegion[numberGoodU_] = pivotMultiplier;
2098 numberInColumn[pivotColumn] = 0;
2099 otherMultiplier = otherMultiplier * pivotMultiplier;
2100 indexRowL[l] = otherRow;
2101 elementL[l] = otherMultiplier;
2102 //take out of row list
2103 CoinBigIndex start = startRowU[otherRow];
2104 CoinBigIndex end = start + numberSave;
2105 CoinBigIndex where = start;
2106 int * indexColumnU = indexColumnU_.array();
2107
2108 while ( indexColumnU[where] != pivotColumn ) {
2109 where++;
2110 } /* endwhile */
2111 assert ( where < end );
2112 end--;
2113 indexColumnU[where] = indexColumnU[end];
2114 int numberAdded = 0;
2115 int numberDeleted = 0;
2116
2117 //pack down and move to work
2118 int j;
2119 const int * nextCount = nextCount_.array();
2120 int * nextColumn = nextColumn_.array();
2121
2122 for ( j = startRow; j < endRow; j++ ) {
2123 int iColumn = indexColumnU[j];
2124
2125 if ( iColumn != pivotColumn ) {
2126 CoinBigIndex startColumn = startColumnU[iColumn];
2127 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
2128 int iRow = indexRowU[startColumn];
2129 CoinFactorizationDouble value = elementU[startColumn];
2130 double largest;
2131 bool foundOther = false;
2132
2133 //leave room for pivot
2134 CoinBigIndex put = startColumn + 1;
2135 CoinBigIndex positionLargest = -1;
2136 CoinFactorizationDouble thisPivotValue = 0.0;
2137 CoinFactorizationDouble otherElement = 0.0;
2138 CoinFactorizationDouble nextValue = elementU[put];
2139 int nextIRow = indexRowU[put];
2140
2141 //compress column and find largest not updated
2142 if ( iRow != pivotRow ) {
2143 if ( iRow != otherRow ) {
2144 largest = fabs ( value );
2145 elementU[put] = value;
2146 indexRowU[put] = iRow;
2147 positionLargest = put;
2148 put++;
2149 CoinBigIndex i;
2150 for ( i = startColumn + 1; i < endColumn; i++ ) {
2151 iRow = nextIRow;
2152 value = nextValue;
2153#ifdef ZEROFAULT
2154 // doesn't matter reading uninitialized but annoys checking
2155 if ( i + 1 < endColumn ) {
2156#endif
2157 nextIRow = indexRowU[i + 1];
2158 nextValue = elementU[i + 1];
2159#ifdef ZEROFAULT
2160 }
2161#endif
2162 if ( iRow != pivotRow ) {
2163 if ( iRow != otherRow ) {
2164 //keep
2165 indexRowU[put] = iRow;
2166 elementU[put] = value;
2167 put++;
2168 } else {
2169 otherElement = value;
2170 foundOther = true;
2171 }
2172 } else {
2173 thisPivotValue = value;
2174 }
2175 }
2176 } else {
2177 otherElement = value;
2178 foundOther = true;
2179 //need to find largest
2180 largest = 0.0;
2181 CoinBigIndex i;
2182 for ( i = startColumn + 1; i < endColumn; i++ ) {
2183 iRow = nextIRow;
2184 value = nextValue;
2185#ifdef ZEROFAULT
2186 // doesn't matter reading uninitialized but annoys checking
2187 if ( i + 1 < endColumn ) {
2188#endif
2189 nextIRow = indexRowU[i + 1];
2190 nextValue = elementU[i + 1];
2191#ifdef ZEROFAULT
2192 }
2193#endif
2194 if ( iRow != pivotRow ) {
2195 //keep
2196 indexRowU[put] = iRow;
2197 elementU[put] = value;
2198 double absValue = fabs ( value );
2199
2200 if ( absValue > largest ) {
2201 largest = absValue;
2202 positionLargest = put;
2203 }
2204 put++;
2205 } else {
2206 thisPivotValue = value;
2207 }
2208 }
2209 }
2210 } else {
2211 //need to find largest
2212 largest = 0.0;
2213 thisPivotValue = value;
2214 CoinBigIndex i;
2215 for ( i = startColumn + 1; i < endColumn; i++ ) {
2216 iRow = nextIRow;
2217 value = nextValue;
2218#ifdef ZEROFAULT
2219 // doesn't matter reading uninitialized but annoys checking
2220 if ( i + 1 < endColumn ) {
2221#endif
2222 nextIRow = indexRowU[i + 1];
2223 nextValue = elementU[i + 1];
2224#ifdef ZEROFAULT
2225 }
2226#endif
2227 if ( iRow != otherRow ) {
2228 //keep
2229 indexRowU[put] = iRow;
2230 elementU[put] = value;
2231 double absValue = fabs ( value );
2232
2233 if ( absValue > largest ) {
2234 largest = absValue;
2235 positionLargest = put;
2236 }
2237 put++;
2238 } else {
2239 otherElement = value;
2240 foundOther = true;
2241 }
2242 }
2243 }
2244 //slot in pivot
2245 elementU[startColumn] = thisPivotValue;
2246 indexRowU[startColumn] = pivotRow;
2247 //clean up counts
2248 startColumn++;
2249 numberInColumn[iColumn] = put - startColumn;
2250 numberInColumnPlus[iColumn]++;
2251 startColumnU[iColumn]++;
2252 otherElement = otherElement - thisPivotValue * otherMultiplier;
2253 double absValue = fabs ( otherElement );
2254
2255 if ( absValue > zeroTolerance_ ) {
2256 if ( !foundOther ) {
2257 //have we space
2258 saveColumn[numberAdded++] = iColumn;
2259 int next = nextColumn[iColumn];
2260 CoinBigIndex space;
2261
2262 space = startColumnU[next] - put - numberInColumnPlus[next];
2263 if ( space <= 0 ) {
2264 //getColumnSpace also moves fixed part
2265 int number = numberInColumn[iColumn];
2266
2267 if ( !getColumnSpace ( iColumn, number + 1 ) ) {
2268 return false;
2269 }
2270 //redo starts
2271 positionLargest =
2272 positionLargest + startColumnU[iColumn] - startColumn;
2273 startColumn = startColumnU[iColumn];
2274 put = startColumn + number;
2275 }
2276 }
2277 elementU[put] = otherElement;
2278 indexRowU[put] = otherRow;
2279 if ( absValue > largest ) {
2280 largest = absValue;
2281 positionLargest = put;
2282 }
2283 put++;
2284 } else {
2285 if ( foundOther ) {
2286 numberDeleted++;
2287 //take out of row list
2288 CoinBigIndex where = start;
2289
2290 while ( indexColumnU[where] != iColumn ) {
2291 where++;
2292 } /* endwhile */
2293 assert ( where < end );
2294 end--;
2295 indexColumnU[where] = indexColumnU[end];
2296 }
2297 }
2298 numberInColumn[iColumn] = put - startColumn;
2299 //move largest
2300 if ( positionLargest >= 0 ) {
2301 value = elementU[positionLargest];
2302 iRow = indexRowU[positionLargest];
2303 elementU[positionLargest] = elementU[startColumn];
2304 indexRowU[positionLargest] = indexRowU[startColumn];
2305 elementU[startColumn] = value;
2306 indexRowU[startColumn] = iRow;
2307 }
2308 //linked list for column
2309 if ( nextCount[iColumn + numberRows_] != -2 ) {
2310 //modify linked list
2311 deleteLink ( iColumn + numberRows_ );
2312 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
2313 }
2314 }
2315 }
2316 //get space for row list
2317 next = nextRow[otherRow];
2318 CoinBigIndex space;
2319
2320 space = startRowU[next] - end;
2321 totalElements_ += numberAdded - numberDeleted;
2322 int number = numberAdded + ( end - start );
2323
2324 if ( space < numberAdded ) {
2325 numberInRow[otherRow] = end - start;
2326 if ( !getRowSpace ( otherRow, number ) ) {
2327 return false;
2328 }
2329 end = startRowU[otherRow] + end - start;
2330 }
2331 // do linked lists and update counts
2332 numberInRow[otherRow] = number;
2333 if ( number != numberSave ) {
2334 deleteLink ( otherRow );
2335 addLink ( otherRow, number );
2336 }
2337 for ( j = 0; j < numberAdded; j++ ) {
2338 indexColumnU[end++] = saveColumn[j];
2339 }
2340 //modify linked list for pivots
2341 deleteLink ( pivotRow );
2342 deleteLink ( pivotColumn + numberRows_ );
2343 return true;
2344}
2345void
2346CoinFactorization::setPersistenceFlag(int flag)
2347{
2348 persistenceFlag_=flag;
2349 workArea_.setPersistence(flag,maximumRowsExtra_+1);
2350 workArea2_.setPersistence(flag,maximumRowsExtra_+1);
2351 pivotColumn_.setPersistence(flag,maximumColumnsExtra_+1);
2352 permute_.setPersistence(flag,maximumRowsExtra_+1);
2353 pivotColumnBack_.setPersistence(flag,maximumRowsExtra_+1);
2354 permuteBack_.setPersistence(flag,maximumRowsExtra_+1);
2355 nextRow_.setPersistence(flag,maximumRowsExtra_+1);
2356 startRowU_.setPersistence(flag,maximumRowsExtra_+1);
2357 numberInRow_.setPersistence(flag,maximumRowsExtra_+1);
2358 numberInColumn_.setPersistence(flag,maximumColumnsExtra_+1);
2359 numberInColumnPlus_.setPersistence(flag,maximumColumnsExtra_+1);
2360 firstCount_.setPersistence(flag,CoinMax(biggerDimension_+2,maximumRowsExtra_+1));
2361 nextCount_.setPersistence(flag,numberRows_+numberColumns_);
2362 lastCount_.setPersistence(flag,numberRows_+numberColumns_);
2363 nextColumn_.setPersistence(flag,maximumColumnsExtra_+1);
2364 lastColumn_.setPersistence(flag,maximumColumnsExtra_+1);
2365 lastRow_.setPersistence(flag,maximumRowsExtra_+1);
2366 markRow_.setPersistence(flag,numberRows_);
2367 saveColumn_.setPersistence(flag,numberColumns_);
2368 indexColumnU_.setPersistence(flag,lengthAreaU_);
2369 pivotRowL_.setPersistence(flag,numberRows_+1);
2370 pivotRegion_.setPersistence(flag,maximumRowsExtra_+1);
2371 elementU_.setPersistence(flag,lengthAreaU_);
2372 indexRowU_.setPersistence(flag,lengthAreaU_);
2373 startColumnU_.setPersistence(flag,maximumColumnsExtra_+1);
2374 convertRowToColumnU_.setPersistence( flag, lengthAreaU_ );
2375 elementL_.setPersistence( flag, lengthAreaL_ );
2376 indexRowL_.setPersistence( flag, lengthAreaL_ );
2377 startColumnL_.setPersistence( flag, numberRows_ + 1 );
2378 startColumnR_.setPersistence( flag, maximumPivots_ + 1 + maximumColumnsExtra_ + 1 );
2379 startRowL_.setPersistence( flag,0 );
2380 indexColumnL_.setPersistence( flag, 0 );
2381 elementByRowL_.setPersistence( flag, 0 );
2382 sparse_.setPersistence( flag, 0 );
2383}
2384// Delete all stuff
2385void
2386CoinFactorization::almostDestructor()
2387{
2388 gutsOfDestructor(1);
2389}
2390