1 /* $Id: CoinSimpFactorization.cpp 1448 2011-06-19 15:34:41Z stefan $ */
2 // Copyright (C) 2008, 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 #include "CoinUtilsConfig.h"
7
8 #include <cassert>
9 #include "CoinPragma.hpp"
10 #include "CoinSimpFactorization.hpp"
11 #include "CoinIndexedVector.hpp"
12 #include "CoinHelperFunctions.hpp"
13 #include "CoinPackedMatrix.hpp"
14 #include "CoinFinite.hpp"
15 #include <stdio.h>
16
17
18
19 #define ARRAY 0
20
21
22
FactorPointers(int numRows,int numColumns,int * UrowLengths_,int * UcolLengths_)23 FactorPointers::FactorPointers( int numRows, int numColumns,
24 int *UrowLengths_,
25 int *UcolLengths_ ){
26
27 rowMax = new double[numRows];
28 double *current=rowMax;
29 const double *end=current+numRows;
30 for ( ; current!=end; ++current ) *current=-1.0;
31
32 firstRowKnonzeros = new int[numRows+1];
33 CoinFillN(firstRowKnonzeros, numRows+1, -1 );
34
35 prevRow = new int[numRows];
36 nextRow = new int[numRows];
37 firstColKnonzeros = new int[numRows+1];
38 memset(firstColKnonzeros, -1, (numRows+1)*sizeof(int) );
39
40 prevColumn = new int[numColumns];
41 nextColumn = new int[numColumns];
42 newCols = new int[numRows];
43
44
45 for ( int i=numRows-1; i>=0; --i ){
46 int length=UrowLengths_[i];
47 prevRow[i]=-1;
48 nextRow[i]=firstRowKnonzeros[length];
49 if ( nextRow[i]!=-1 ) prevRow[nextRow[i]]=i;
50 firstRowKnonzeros[length]=i;
51 }
52 for ( int i=numColumns-1; i>=0; --i ){
53 int length=UcolLengths_[i];
54 prevColumn[i]=-1;
55 nextColumn[i]=firstColKnonzeros[length];
56 if ( nextColumn[i]!=-1 ) prevColumn[nextColumn[i]]=i;
57 firstColKnonzeros[length]=i;
58 }
59 }
~FactorPointers()60 FactorPointers::~ FactorPointers(){
61 delete [] rowMax;
62 delete [] firstRowKnonzeros;
63 delete [] prevRow;
64 delete [] nextRow;
65 delete [] firstColKnonzeros;
66 delete [] prevColumn;
67 delete [] nextColumn;
68 delete [] newCols;
69 }
70
71
72 //:class CoinSimpFactorization. Deals with Factorization and Updates
73 // CoinSimpFactorization. Constructor
CoinSimpFactorization()74 CoinSimpFactorization::CoinSimpFactorization ( )
75 : CoinOtherFactorization()
76 {
77 gutsOfInitialize();
78 }
79
80 /// Copy constructor
CoinSimpFactorization(const CoinSimpFactorization & other)81 CoinSimpFactorization::CoinSimpFactorization ( const CoinSimpFactorization &other)
82 : CoinOtherFactorization(other)
83 {
84 gutsOfInitialize();
85 gutsOfCopy(other);
86 }
87 // Clone
88 CoinOtherFactorization *
clone() const89 CoinSimpFactorization::clone() const
90 {
91 return new CoinSimpFactorization(*this);
92 }
93 /// The real work of constructors etc
gutsOfDestructor()94 void CoinSimpFactorization::gutsOfDestructor()
95 {
96 delete [] elements_;
97 delete [] pivotRow_;
98 delete [] workArea_;
99 elements_ = NULL;
100 pivotRow_ = NULL;
101 workArea_ = NULL;
102 numberRows_ = 0;
103 numberColumns_ = 0;
104 numberGoodU_ = 0;
105 status_ = -1;
106 maximumRows_=0;
107 maximumSpace_=0;
108 numberSlacks_=0;
109 firstNumberSlacks_=0;
110 //
111 delete [] denseVector_;
112 delete [] workArea2_;
113 delete [] workArea3_;
114 delete [] vecLabels_;
115 delete [] indVector_;
116
117 delete [] auxVector_;
118 delete [] auxInd_;
119
120 delete [] vecKeep_;
121 delete [] indKeep_;
122
123 delete [] LrowStarts_;
124 delete [] LrowLengths_;
125 delete [] Lrows_;
126 delete [] LrowInd_;
127
128
129 delete [] LcolStarts_;
130 delete [] LcolLengths_;
131 delete [] Lcolumns_;
132 delete [] LcolInd_;
133
134 delete [] UrowStarts_;
135 delete [] UrowLengths_;
136 #ifdef COIN_SIMP_CAPACITY
137 delete [] UrowCapacities_;
138 #endif
139 delete [] Urows_;
140 delete [] UrowInd_;
141
142 delete [] prevRowInU_;
143 delete [] nextRowInU_;
144
145 delete [] UcolStarts_;
146 delete [] UcolLengths_;
147 #ifdef COIN_SIMP_CAPACITY
148 delete [] UcolCapacities_;
149 #endif
150 delete [] Ucolumns_;
151 delete [] UcolInd_;
152 delete [] prevColInU_;
153 delete [] nextColInU_;
154 delete [] colSlack_;
155
156 delete [] invOfPivots_;
157
158 delete [] colOfU_;
159 delete [] colPosition_;
160 delete [] rowOfU_;
161 delete [] rowPosition_;
162 delete [] secRowOfU_;
163 delete [] secRowPosition_;
164
165 delete [] EtaPosition_;
166 delete [] EtaStarts_;
167 delete [] EtaLengths_;
168 delete [] EtaInd_;
169 delete [] Eta_;
170
171 denseVector_=NULL;
172 workArea2_=NULL;
173 workArea3_=NULL;
174 vecLabels_=NULL;
175 indVector_=NULL;
176
177 auxVector_=NULL;
178 auxInd_=NULL;
179
180 vecKeep_=NULL;
181 indKeep_=NULL;
182
183 LrowStarts_=NULL;
184 LrowLengths_=NULL;
185 Lrows_=NULL;
186 LrowInd_=NULL;
187
188
189 LcolStarts_=NULL;
190 LcolLengths_=NULL;
191 Lcolumns_=NULL;
192 LcolInd_=NULL;
193
194 UrowStarts_=NULL;
195 UrowLengths_=NULL;
196 #ifdef COIN_SIMP_CAPACITY
197 UrowCapacities_=NULL;
198 #endif
199 Urows_=NULL;
200 UrowInd_=NULL;
201
202 prevRowInU_=NULL;
203 nextRowInU_=NULL;
204
205 UcolStarts_=NULL;
206 UcolLengths_=NULL;
207 #ifdef COIN_SIMP_CAPACITY
208 UcolCapacities_=NULL;
209 #endif
210 Ucolumns_=NULL;
211 UcolInd_=NULL;
212 prevColInU_=NULL;
213 nextColInU_=NULL;
214 colSlack_=NULL;
215
216 invOfPivots_=NULL;
217
218 colOfU_=NULL;
219 colPosition_=NULL;
220 rowOfU_=NULL;
221 rowPosition_=NULL;
222 secRowOfU_=NULL;
223 secRowPosition_=NULL;
224
225 EtaPosition_=NULL;
226 EtaStarts_=NULL;
227 EtaLengths_=NULL;
228 EtaInd_=NULL;
229 Eta_=NULL;
230 }
gutsOfInitialize()231 void CoinSimpFactorization::gutsOfInitialize()
232 {
233 pivotTolerance_ = 1.0e-1;
234 zeroTolerance_ = 1.0e-13;
235 #ifndef COIN_FAST_CODE
236 slackValue_ = -1.0;
237 #endif
238 maximumPivots_=200;
239 relaxCheck_=1.0;
240 numberRows_ = 0;
241 numberColumns_ = 0;
242 numberGoodU_ = 0;
243 status_ = -1;
244 numberPivots_ = 0;
245 maximumRows_=0;
246 maximumSpace_=0;
247 numberSlacks_=0;
248 firstNumberSlacks_=0;
249 elements_ = NULL;
250 pivotRow_ = NULL;
251 workArea_ = NULL;
252
253 denseVector_=NULL;
254 workArea2_=NULL;
255 workArea3_=NULL;
256 vecLabels_=NULL;
257 indVector_=NULL;
258
259 auxVector_=NULL;
260 auxInd_=NULL;
261
262 vecKeep_=NULL;
263 indKeep_=NULL;
264
265 LrowStarts_=NULL;
266 LrowLengths_=NULL;
267 Lrows_=NULL;
268 LrowInd_=NULL;
269
270
271 LcolStarts_=NULL;
272 LcolLengths_=NULL;
273 Lcolumns_=NULL;
274 LcolInd_=NULL;
275
276 UrowStarts_=NULL;
277 UrowLengths_=NULL;
278 #ifdef COIN_SIMP_CAPACITY
279 UrowCapacities_=NULL;
280 #endif
281 Urows_=NULL;
282 UrowInd_=NULL;
283
284 prevRowInU_=NULL;
285 nextRowInU_=NULL;
286
287 UcolStarts_=NULL;
288 UcolLengths_=NULL;
289 #ifdef COIN_SIMP_CAPACITY
290 UcolCapacities_=NULL;
291 #endif
292 Ucolumns_=NULL;
293 UcolInd_=NULL;
294 prevColInU_=NULL;
295 nextColInU_=NULL;
296 colSlack_=NULL;
297
298 invOfPivots_=NULL;
299
300 colOfU_=NULL;
301 colPosition_=NULL;
302 rowOfU_=NULL;
303 rowPosition_=NULL;
304 secRowOfU_=NULL;
305 secRowPosition_=NULL;
306
307 EtaPosition_=NULL;
308 EtaStarts_=NULL;
309 EtaLengths_=NULL;
310 EtaInd_=NULL;
311 Eta_=NULL;
312 }
initialSomeNumbers()313 void CoinSimpFactorization::initialSomeNumbers(){
314 keepSize_=-1;
315 LrowSize_=-1;
316 // LrowCap_ in allocateSomeArrays
317 LcolSize_=-1;
318 // LcolCap_ in allocateSomeArrays
319 // UrowMaxCap_ in allocateSomeArrays
320 UrowEnd_=-1;
321 firstRowInU_=-1;
322 lastRowInU_=-1;
323 firstColInU_=-1;
324 lastColInU_=-1;
325 //UcolMaxCap_ in allocateSomeArrays
326 UcolEnd_=-1;
327
328
329 EtaSize_=0;
330 lastEtaRow_=-1;
331 //maxEtaRows_ in allocateSomeArrays
332 //EtaMaxCap_ in allocateSomeArrays
333
334 // minIncrease_ in allocateSomeArrays
335 updateTol_=1.0e12;
336
337 doSuhlHeuristic_=true;
338 maxU_=-1.0;
339 maxGrowth_=1.e12;
340 maxA_=-1.0;
341 pivotCandLimit_=4;
342 minIncrease_=10;
343
344 }
345
346 // ~CoinSimpFactorization. Destructor
~CoinSimpFactorization()347 CoinSimpFactorization::~CoinSimpFactorization ( )
348 {
349 gutsOfDestructor();
350 }
351 // =
operator =(const CoinSimpFactorization & other)352 CoinSimpFactorization & CoinSimpFactorization::operator = ( const CoinSimpFactorization & other ) {
353 if (this != &other) {
354 gutsOfDestructor();
355 gutsOfInitialize();
356 gutsOfCopy(other);
357 }
358 return *this;
359 }
gutsOfCopy(const CoinSimpFactorization & other)360 void CoinSimpFactorization::gutsOfCopy(const CoinSimpFactorization &other)
361 {
362 pivotTolerance_ = other.pivotTolerance_;
363 zeroTolerance_ = other.zeroTolerance_;
364 #ifndef COIN_FAST_CODE
365 slackValue_ = other.slackValue_;
366 #endif
367 relaxCheck_ = other.relaxCheck_;
368 numberRows_ = other.numberRows_;
369 numberColumns_ = other.numberColumns_;
370 maximumRows_ = other.maximumRows_;
371 maximumSpace_ = other.maximumSpace_;
372 numberGoodU_ = other.numberGoodU_;
373 maximumPivots_ = other.maximumPivots_;
374 numberPivots_ = other.numberPivots_;
375 factorElements_ = other.factorElements_;
376 status_ = other.status_;
377 numberSlacks_ = other.numberSlacks_;
378 firstNumberSlacks_ = other.firstNumberSlacks_;
379 if (other.pivotRow_) {
380 pivotRow_ = new int [2*maximumRows_+maximumPivots_];
381 memcpy(pivotRow_,other.pivotRow_,(2*maximumRows_+numberPivots_)*sizeof(int));
382 elements_ = new CoinFactorizationDouble [maximumSpace_];
383 memcpy(elements_,other.elements_,(maximumRows_+numberPivots_)*maximumRows_*sizeof(CoinFactorizationDouble));
384 workArea_ = new CoinFactorizationDouble [maximumRows_];
385 } else {
386 elements_ = NULL;
387 pivotRow_ = NULL;
388 workArea_ = NULL;
389 }
390
391 keepSize_ = other.keepSize_;
392
393 LrowSize_ = other.LrowSize_;
394 LrowCap_ = other.LrowCap_;
395
396 LcolSize_ = other.LcolSize_;
397 LcolCap_ = other.LcolCap_;
398
399 UrowMaxCap_ = other.UrowMaxCap_;
400 UrowEnd_ = other.UrowEnd_;
401 firstRowInU_ = other.firstRowInU_;
402 lastRowInU_ = other.lastRowInU_;
403
404 firstColInU_ = other.firstColInU_;
405 lastColInU_ = other.lastColInU_;
406 UcolMaxCap_ = other.UcolMaxCap_;
407 UcolEnd_ = other.UcolEnd_;
408
409 EtaSize_ = other.EtaSize_;
410 lastEtaRow_ = other.lastEtaRow_;
411 maxEtaRows_ = other.maxEtaRows_;
412 EtaMaxCap_ = other.EtaMaxCap_;
413
414
415 minIncrease_ = other.minIncrease_;
416 updateTol_ = other.updateTol_;
417
418
419
420 if (other.denseVector_)
421 {
422 denseVector_ = new double[maximumRows_];
423 memcpy(denseVector_,other.denseVector_,maximumRows_*sizeof(double));
424 } else denseVector_=NULL;
425 if (other.workArea2_)
426 {
427 workArea2_ = new double[maximumRows_];
428 memcpy(workArea2_,other.workArea2_,maximumRows_*sizeof(double));
429 } else workArea2_=NULL;
430 if (other.workArea3_)
431 {
432 workArea3_ = new double[maximumRows_];
433 memcpy(workArea3_,other.workArea3_,maximumRows_*sizeof(double));
434 } else workArea3_=NULL;
435 if (other.vecLabels_)
436 {
437 vecLabels_ = new int[maximumRows_];
438 memcpy(vecLabels_,other.vecLabels_,maximumRows_*sizeof(int));
439 } else vecLabels_=NULL;
440 if (other.indVector_)
441 {
442 indVector_ = new int[maximumRows_];
443 memcpy(indVector_ ,other.indVector_ , maximumRows_ *sizeof(int ));
444 } else indVector_=NULL;
445
446 if (other.auxVector_)
447 {
448 auxVector_ = new double[maximumRows_];
449 memcpy(auxVector_ ,other.auxVector_ , maximumRows_ *sizeof(double ));
450 } else auxVector_=NULL;
451 if (other.auxInd_)
452 {
453 auxInd_ = new int[maximumRows_];
454 memcpy(auxInd_ , other.auxInd_, maximumRows_ *sizeof(int));
455 } else auxInd_=NULL;
456
457 if (other.vecKeep_)
458 {
459 vecKeep_ = new double[maximumRows_];
460 memcpy(vecKeep_ ,other.vecKeep_ , maximumRows_ *sizeof(double));
461 } else vecKeep_=NULL;
462 if (other.indKeep_)
463 {
464 indKeep_ = new int[maximumRows_];
465 memcpy(indKeep_ , other.indKeep_, maximumRows_ *sizeof(int));
466 } else indKeep_=NULL;
467 if (other.LrowStarts_)
468 {
469 LrowStarts_ = new int[maximumRows_];
470 memcpy(LrowStarts_ , other.LrowStarts_, maximumRows_ *sizeof(int));
471 } else LrowStarts_=NULL;
472 if (other.LrowLengths_)
473 {
474 LrowLengths_ = new int[maximumRows_];
475 memcpy(LrowLengths_ , other.LrowLengths_ , maximumRows_ *sizeof(int));
476 } else LrowLengths_=NULL;
477 if (other.Lrows_)
478 {
479 Lrows_ = new double[other.LrowCap_];
480 memcpy(Lrows_ , other.Lrows_, other.LrowCap_*sizeof(double));
481 } else Lrows_=NULL;
482 if (other.LrowInd_)
483 {
484 LrowInd_ = new int[other.LrowCap_];
485 memcpy(LrowInd_, other.LrowInd_, other.LrowCap_*sizeof(int));
486 } else LrowInd_=NULL;
487
488 if (other.LcolStarts_)
489 {
490 LcolStarts_ = new int[maximumRows_];
491 memcpy(LcolStarts_ ,other.LcolStarts_ , maximumRows_*sizeof(int));
492 } else LcolStarts_=NULL;
493 if (other.LcolLengths_)
494 {
495 LcolLengths_ = new int[maximumRows_];
496 memcpy(LcolLengths_ , other.LcolLengths_ , maximumRows_ *sizeof(int));
497 } else LcolLengths_=NULL;
498 if (other.Lcolumns_)
499 {
500 Lcolumns_ = new double[other.LcolCap_];
501 memcpy(Lcolumns_ ,other.Lcolumns_ , other.LcolCap_ *sizeof(double));
502 } else Lcolumns_=NULL;
503 if (other.LcolInd_)
504 {
505 LcolInd_ = new int[other.LcolCap_];
506 memcpy(LcolInd_ , other.LcolInd_, other.LcolCap_*sizeof(int));
507 } else LcolInd_=NULL;
508
509 if (other.UrowStarts_)
510 {
511 UrowStarts_ = new int[maximumRows_];
512 memcpy(UrowStarts_ ,other.UrowStarts_ , maximumRows_ *sizeof(int));
513 } else UrowStarts_=NULL;
514 if (other.UrowLengths_)
515 {
516 UrowLengths_ = new int[maximumRows_];
517 memcpy(UrowLengths_ ,other.UrowLengths_ , maximumRows_ *sizeof(int));
518 } else UrowLengths_=NULL;
519 #ifdef COIN_SIMP_CAPACITY
520 if (other.UrowCapacities_)
521 {
522 UrowCapacities_ = new int[maximumRows_];
523 memcpy(UrowCapacities_ ,other.UrowCapacities_ , maximumRows_ *sizeof(int));
524 } else UrowCapacities_=NULL;
525 #endif
526 if (other.Urows_)
527 {
528 Urows_ = new double[other.UrowMaxCap_];
529 memcpy(Urows_ ,other.Urows_ , other.UrowMaxCap_ *sizeof(double));
530 } else Urows_=NULL;
531 if (other.UrowInd_)
532 {
533 UrowInd_ = new int[other.UrowMaxCap_];
534 memcpy(UrowInd_,other.UrowInd_, other.UrowMaxCap_*sizeof(int));
535 } else UrowInd_=NULL;
536
537 if (other.prevRowInU_)
538 {
539 prevRowInU_ = new int[maximumRows_];
540 memcpy(prevRowInU_ , other.prevRowInU_, maximumRows_*sizeof(int));
541 } else prevRowInU_=NULL;
542 if (other.nextRowInU_)
543 {
544 nextRowInU_ = new int[maximumRows_];
545 memcpy(nextRowInU_, other.nextRowInU_, maximumRows_*sizeof(int));
546 } else nextRowInU_=NULL;
547
548 if (other.UcolStarts_)
549 {
550 UcolStarts_ = new int[maximumRows_];
551 memcpy(UcolStarts_ , other.UcolStarts_, maximumRows_*sizeof(int));
552 } else UcolStarts_=NULL;
553 if (other.UcolLengths_)
554 {
555 UcolLengths_ = new int[maximumRows_];
556 memcpy(UcolLengths_ , other.UcolLengths_, maximumRows_*sizeof(int));
557 } else UcolLengths_=NULL;
558 #ifdef COIN_SIMP_CAPACITY
559 if (other.UcolCapacities_)
560 {
561 UcolCapacities_ = new int[maximumRows_];
562 memcpy(UcolCapacities_ ,other.UcolCapacities_ , maximumRows_*sizeof(int));
563 } else UcolCapacities_=NULL;
564 #endif
565 if (other.Ucolumns_)
566 {
567 Ucolumns_ = new double[other.UcolMaxCap_];
568 memcpy(Ucolumns_ ,other.Ucolumns_ , other.UcolMaxCap_*sizeof(double));
569 } else Ucolumns_=NULL;
570 if (other.UcolInd_)
571 {
572 UcolInd_ = new int[other.UcolMaxCap_];
573 memcpy(UcolInd_ , other.UcolInd_ , other.UcolMaxCap_*sizeof(int));
574 } else UcolInd_=NULL;
575 if (other.prevColInU_)
576 {
577 prevColInU_ = new int[maximumRows_];
578 memcpy(prevColInU_ , other.prevColInU_ , maximumRows_*sizeof(int));
579 } else prevColInU_=NULL;
580 if (other.nextColInU_)
581 {
582 nextColInU_ = new int[maximumRows_];
583 memcpy(nextColInU_ ,other.nextColInU_ , maximumRows_*sizeof(int));
584 } else nextColInU_=NULL;
585 if (other.colSlack_)
586 {
587 colSlack_ = new int[maximumRows_];
588 memcpy(colSlack_, other.colSlack_, maximumRows_*sizeof(int));
589 }
590
591
592 if (other.invOfPivots_)
593 {
594 invOfPivots_ = new double[maximumRows_];
595 memcpy(invOfPivots_ , other.invOfPivots_, maximumRows_*sizeof(double));
596 } else invOfPivots_=NULL;
597
598 if (other.colOfU_)
599 {
600 colOfU_ = new int[maximumRows_];
601 memcpy(colOfU_ , other.colOfU_, maximumRows_*sizeof(int));
602 } else colOfU_=NULL;
603 if (other.colPosition_)
604 {
605 colPosition_ = new int[maximumRows_];
606 memcpy(colPosition_, other.colPosition_, maximumRows_*sizeof(int));
607 } else colPosition_=NULL;
608 if (other.rowOfU_)
609 {
610 rowOfU_ = new int[maximumRows_];
611 memcpy(rowOfU_ , other.rowOfU_, maximumRows_*sizeof(int));
612 } else rowOfU_=NULL;
613 if (other.rowPosition_)
614 {
615 rowPosition_ = new int[maximumRows_];
616 memcpy(rowPosition_ , other.rowPosition_, maximumRows_*sizeof(int));
617 } else rowPosition_=NULL;
618 if (other.secRowOfU_)
619 {
620 secRowOfU_ = new int[maximumRows_];
621 memcpy(secRowOfU_ , other.secRowOfU_, maximumRows_*sizeof(int));
622 } else secRowOfU_=NULL;
623 if (other.secRowPosition_)
624 {
625 secRowPosition_ = new int[maximumRows_];
626 memcpy(secRowPosition_ , other.secRowPosition_, maximumRows_*sizeof(int));
627 } else secRowPosition_=NULL;
628
629 if (other.EtaPosition_)
630 {
631 EtaPosition_ = new int[other.maxEtaRows_];
632 memcpy(EtaPosition_ ,other.EtaPosition_ , other.maxEtaRows_ *sizeof(int));
633 } else EtaPosition_=NULL;
634 if (other.EtaStarts_)
635 {
636 EtaStarts_ = new int[other.maxEtaRows_];
637 memcpy(EtaStarts_, other.EtaStarts_, other.maxEtaRows_*sizeof(int));
638 } else EtaStarts_=NULL;
639 if (other.EtaLengths_)
640 {
641 EtaLengths_ = new int[other.maxEtaRows_];
642 memcpy(EtaLengths_, other.EtaLengths_, other.maxEtaRows_*sizeof(int));
643 } else EtaLengths_=NULL;
644 if (other.EtaInd_)
645 {
646 EtaInd_ = new int[other.EtaMaxCap_];
647 memcpy(EtaInd_, other.EtaInd_, other.EtaMaxCap_*sizeof(int));
648 } else EtaInd_=NULL;
649 if (other.Eta_)
650 {
651 Eta_ = new double[other.EtaMaxCap_];
652 memcpy(Eta_ , other.Eta_, other.EtaMaxCap_*sizeof(double));
653 } else Eta_=NULL;
654
655
656
657 doSuhlHeuristic_ = other.doSuhlHeuristic_;
658 maxU_ = other.maxU_;
659 maxGrowth_ = other.maxGrowth_;
660 maxA_ = other.maxA_;
661 pivotCandLimit_ = other.pivotCandLimit_;
662 }
663
664 // getAreas. Gets space for a factorization
665 //called by constructors
666 void
getAreas(int numberOfRows,int numberOfColumns,CoinBigIndex,CoinBigIndex)667 CoinSimpFactorization::getAreas ( int numberOfRows,
668 int numberOfColumns,
669 CoinBigIndex ,
670 CoinBigIndex )
671 {
672
673 numberRows_ = numberOfRows;
674 numberColumns_ = numberOfColumns;
675 CoinBigIndex size = numberRows_*(numberRows_+CoinMax(maximumPivots_,(numberRows_+1)>>1));
676 if (size>maximumSpace_) {
677 delete [] elements_;
678 elements_ = new CoinFactorizationDouble [size];
679 maximumSpace_ = size;
680 }
681 if (numberRows_>maximumRows_) {
682 maximumRows_ = numberRows_;
683 delete [] pivotRow_;
684 delete [] workArea_;
685 pivotRow_ = new int [2*maximumRows_+maximumPivots_];
686 workArea_ = new CoinFactorizationDouble [maximumRows_];
687 allocateSomeArrays();
688 }
689 }
690
691 // preProcess.
692 void
preProcess()693 CoinSimpFactorization::preProcess ()
694 {
695 CoinBigIndex put = numberRows_*numberRows_;
696 int *indexRow = reinterpret_cast<int *> (elements_+put);
697 CoinBigIndex * starts = reinterpret_cast<CoinBigIndex *> (pivotRow_);
698 initialSomeNumbers();
699
700 // compute sizes for Urows_ and Ucolumns_
701 //for ( int row=0; row < numberRows_; ++row )
702 //UrowLengths_[row]=0;
703 int k=0;
704 for ( int column=0; column < numberColumns_; ++column ){
705 UcolStarts_[column]=k;
706 //for (CoinBigIndex j=starts[column];j<starts[column+1];j++) {
707 // int iRow = indexRow[j];
708 // ++UrowLengths_[iRow];
709 // }
710 UcolLengths_[column]=starts[column+1]-starts[column];
711 #ifdef COIN_SIMP_CAPACITY
712 UcolCapacities_[column]=numberRows_;
713 #endif
714 //k+=UcolLengths_[column]+minIncrease_;
715 k+=numberRows_;
716 }
717
718 // create space for Urows_
719 k=0;
720 for ( int row=0; row < numberRows_; ++row ){
721 prevRowInU_[row]=row-1;
722 nextRowInU_[row]=row+1;
723 UrowStarts_[row]=k;
724 //k+=UrowLengths_[row]+minIncrease_;
725 k+=numberRows_;
726 UrowLengths_[row]=0;
727 #ifdef COIN_SIMP_CAPACITY
728 UrowCapacities_[row]=numberRows_;
729 #endif
730 }
731 UrowEnd_=k;
732 nextRowInU_[numberRows_-1]=-1;
733 firstRowInU_=0;
734 lastRowInU_=numberRows_-1;
735 maxA_=-1.0;
736 // build Ucolumns_ and Urows_
737 int colBeg, colEnd;
738 for ( int column=0; column < numberColumns_; ++column ){
739 prevColInU_[column]=column-1;
740 nextColInU_[column]=column+1;
741 k=0;
742 colBeg=starts[column];
743 colEnd=starts[column+1];
744 // identify slacks
745 if ( colEnd == colBeg+1 && elements_[colBeg]==slackValue_ )
746 colSlack_[column]=1;
747 else colSlack_[column]=0;
748 //
749 for (int j=colBeg; j < colEnd; j++) {
750 // Ucolumns_
751 int iRow = indexRow[j];
752 UcolInd_[UcolStarts_[column] + k]=iRow;
753 ++k;
754 // Urow_
755 int ind=UrowStarts_[iRow]+UrowLengths_[iRow];
756 UrowInd_[ind]=column;
757 Urows_[ind]=elements_[j];
758 //maxA_=CoinMax( maxA_, fabs(Urows_[ind]) );
759 ++UrowLengths_[iRow];
760 }
761 }
762 nextColInU_[numberColumns_-1]=-1;
763 firstColInU_=0;
764 lastColInU_=numberColumns_-1;
765
766 // Initialize L
767 LcolSize_=0;
768 //LcolCap_=numberRows_*minIncrease_;
769 memset(LrowStarts_,-1,numberRows_ * sizeof(int));
770 memset(LrowLengths_,0,numberRows_ * sizeof(int));
771 memset(LcolStarts_,-1,numberRows_ * sizeof(int));
772 memset(LcolLengths_,0,numberRows_ * sizeof(int));
773
774 // initialize permutations
775 for ( int i=0; i<numberRows_; ++i ){
776 rowOfU_[i]=i;
777 rowPosition_[i]=i;
778 }
779 for ( int i=0; i<numberColumns_; ++i ){
780 colOfU_[i]=i;
781 colPosition_[i]=i;
782 }
783 //
784
785 doSuhlHeuristic_=true;
786
787 }
788
factorize(int numberOfRows,int numberOfColumns,const int colStarts[],const int indicesRow[],const double elements[])789 void CoinSimpFactorization::factorize(int numberOfRows,
790 int numberOfColumns,
791 const int colStarts[],
792 const int indicesRow[],
793 const double elements[])
794 {
795 CoinBigIndex maximumL=0;
796 CoinBigIndex maximumU=0;
797 getAreas ( numberOfRows, numberOfColumns, maximumL, maximumU );
798
799 CoinBigIndex put = numberRows_*numberRows_;
800 int *indexRow = reinterpret_cast<int *> (elements_+put);
801 CoinBigIndex * starts = reinterpret_cast<CoinBigIndex *> (pivotRow_);
802 for ( int column=0; column <= numberColumns_; ++column ){
803 starts[column]=colStarts[column];
804 }
805 const int limit=colStarts[numberColumns_];
806 for ( int i=0; i<limit; ++i ){
807 indexRow[i]=indicesRow[i];
808 elements_[i]=elements[i];
809 }
810
811 preProcess();
812 factor();
813 }
814
815 //Does factorization
816 int
factor()817 CoinSimpFactorization::factor ( )
818 {
819 numberPivots_=0;
820 status_= 0;
821
822 FactorPointers pointers(numberRows_, numberColumns_, UrowLengths_, UcolLengths_);
823 int rc=mainLoopFactor (pointers);
824 // rc=0 success
825 if ( rc != 0 ) { // failure
826 status_ = -1;
827 //return status_; // failure
828 }
829 //if ( rc == -3 ) { // numerical instability
830 // status_ = -1;
831 // return status_;
832 //}
833
834 //copyLbyRows();
835 copyUbyColumns();
836 copyRowPermutations();
837 firstNumberSlacks_=numberSlacks_;
838 // row permutations
839 if ( status_==-1 || numberColumns_ < numberRows_ ){
840 for (int j=0;j<numberRows_;j++)
841 pivotRow_[j+numberRows_]=rowOfU_[j];
842 for (int j=0;j<numberRows_;j++) {
843 int k = pivotRow_[j+numberRows_];
844 pivotRow_[k]=j;
845 }
846 }
847 else // no permutations
848 for (int j=0;j<numberRows_;j++){
849 pivotRow_[j]=j;
850 pivotRow_[j+numberRows_]=j;
851 }
852
853 return status_;
854 }
855 //
856
857
858
859 // Makes a non-singular basis by replacing variables
860 void
makeNonSingular(int * sequence,int numberColumns)861 CoinSimpFactorization::makeNonSingular(int * sequence, int numberColumns)
862 {
863 // Replace bad ones by correct slack
864 int * workArea= reinterpret_cast<int *> (workArea_);
865 int i;
866 for ( i=0;i<numberRows_;i++)
867 workArea[i]=-1;
868 for ( i=0;i<numberGoodU_;i++) {
869 int iOriginal = pivotRow_[i+numberRows_];
870 workArea[iOriginal]=i;
871 //workArea[i]=iOriginal;
872 }
873 int lastRow=-1;
874 for ( i=0;i<numberRows_;i++) {
875 if (workArea[i]==-1) {
876 lastRow=i;
877 break;
878 }
879 }
880 assert (lastRow>=0);
881 for ( i=numberGoodU_;i<numberRows_;i++) {
882 assert (lastRow<numberRows_);
883 // Put slack in basis
884 sequence[i]=lastRow+numberColumns;
885 lastRow++;
886 for (;lastRow<numberRows_;lastRow++) {
887 if (workArea[lastRow]==-1)
888 break;
889 }
890 }
891 }
892 // Does post processing on valid factorization - putting variables on correct rows
893 void
postProcess(const int * sequence,int * pivotVariable)894 CoinSimpFactorization::postProcess(const int * sequence, int * pivotVariable)
895 {
896 for (int i=0;i<numberRows_;i++) {
897 int k = sequence[i];
898 pivotVariable[pivotRow_[i+numberRows_]]=k;
899 }
900 }
901 /* Replaces one Column to basis,
902 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
903 If checkBeforeModifying is true will do all accuracy checks
904 before modifying factorization. Whether to set this depends on
905 speed considerations. You could just do this on first iteration
906 after factorization and thereafter re-factorize
907 partial update already in U */
908 int
replaceColumn(CoinIndexedVector *,int pivotRow,double pivotCheck,bool,double)909 CoinSimpFactorization::replaceColumn ( CoinIndexedVector * ,
910 int pivotRow,
911 double pivotCheck ,
912 bool ,
913 double )
914 {
915 if (numberPivots_==maximumPivots_)
916 return 3;
917
918 double pivotValue = pivotCheck;
919 if (fabs(pivotValue)<zeroTolerance_)
920 return 2;
921 int realPivotRow = pivotRow_[pivotRow];
922
923 LUupdate(pivotRow);
924
925 pivotRow_[2*numberRows_+numberPivots_]=realPivotRow;
926 numberPivots_++;
927
928 return 0;
929 }
930 int
updateColumn(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2,bool noPermute) const931 CoinSimpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
932 CoinIndexedVector * regionSparse2,
933 bool noPermute) const
934 {
935 return upColumn(regionSparse, regionSparse2, noPermute, false);
936 }
937 int
updateColumnFT(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2,bool noPermute)938 CoinSimpFactorization::updateColumnFT( CoinIndexedVector * regionSparse,
939 CoinIndexedVector * regionSparse2,
940 bool noPermute)
941 {
942
943 int rc=upColumn(regionSparse, regionSparse2, noPermute, true);
944 return rc;
945 }
946
947
948 int
upColumn(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2,bool,bool save) const949 CoinSimpFactorization::upColumn( CoinIndexedVector * regionSparse,
950 CoinIndexedVector * regionSparse2,
951 bool , bool save) const
952 {
953 assert (numberRows_==numberColumns_);
954 double *region2 = regionSparse2->denseVector ( );
955 int *regionIndex = regionSparse2->getIndices ( );
956 int numberNonZero = regionSparse2->getNumElements ( );
957 double *region=regionSparse->denseVector ( );
958 if (!regionSparse2->packedMode()) {
959 region=regionSparse2->denseVector ( );
960 }
961 else { // packed mode
962 for (int j=0;j<numberNonZero;j++) {
963 region[regionIndex[j]]=region2[j];
964 region2[j]=0.0;
965 }
966 }
967
968 double *solution=workArea2_;
969 ftran(region, solution, save);
970
971 // get nonzeros
972 numberNonZero=0;
973 if (!regionSparse2->packedMode()) {
974 for (int i=0;i<numberRows_;i++) {
975 const double value=solution[i];
976 if ( fabs(value) > zeroTolerance_ ){
977 region[i]=value;
978 regionIndex[numberNonZero++]=i;
979 }
980 else
981 region[i]=0.0;
982 }
983 }
984 else { // packed mode
985 memset(region,0,numberRows_*sizeof(double));
986 for (int i=0;i<numberRows_;i++) {
987 const double value=solution[i];
988 if ( fabs(value) > zeroTolerance_ ){
989 region2[numberNonZero] = value;
990 regionIndex[numberNonZero++]=i;
991 }
992 }
993 }
994 regionSparse2->setNumElements(numberNonZero);
995 return 0;
996 }
997
998
999 int
updateTwoColumnsFT(CoinIndexedVector * regionSparse1,CoinIndexedVector * regionSparse2,CoinIndexedVector * regionSparse3,bool)1000 CoinSimpFactorization::updateTwoColumnsFT(CoinIndexedVector * regionSparse1,
1001 CoinIndexedVector * regionSparse2,
1002 CoinIndexedVector * regionSparse3,
1003 bool )
1004 {
1005 assert (numberRows_==numberColumns_);
1006
1007 double *region2 = regionSparse2->denseVector ( );
1008 int *regionIndex2 = regionSparse2->getIndices ( );
1009 int numberNonZero2 = regionSparse2->getNumElements ( );
1010
1011 double *vec1=regionSparse1->denseVector ( );
1012 if (!regionSparse2->packedMode()) {
1013 vec1=regionSparse2->denseVector ( );
1014 }
1015 else { // packed mode
1016 for (int j=0;j<numberNonZero2;j++) {
1017 vec1[regionIndex2[j]]=region2[j];
1018 region2[j]=0.0;
1019 }
1020 }
1021 //
1022 double *region3 = regionSparse3->denseVector ( );
1023 int *regionIndex3 = regionSparse3->getIndices ( );
1024 int numberNonZero3 = regionSparse3->getNumElements ( );
1025 double *vec2=auxVector_;
1026 if (!regionSparse3->packedMode()) {
1027 vec2=regionSparse3->denseVector ( );
1028 }
1029 else { // packed mode
1030 memset(vec2,0,numberRows_*sizeof(double));
1031 for (int j=0;j<numberNonZero3;j++) {
1032 vec2[regionIndex3[j]]=region3[j];
1033 region3[j]=0.0;
1034 }
1035 }
1036
1037 double *solution1=workArea2_;
1038 double *solution2=workArea3_;
1039 ftran2(vec1, solution1, vec2, solution2);
1040
1041 // get nonzeros
1042 numberNonZero2=0;
1043 if (!regionSparse2->packedMode()) {
1044 double value;
1045 for (int i=0;i<numberRows_;i++) {
1046 value=solution1[i];
1047 if ( fabs(value) > zeroTolerance_ ){
1048 vec1[i]=value;
1049 regionIndex2[numberNonZero2++]=i;
1050 }
1051 else
1052 vec1[i]=0.0;
1053 }
1054 }
1055 else { // packed mode
1056 double value;
1057 for (int i=0;i<numberRows_;i++) {
1058 vec1[i]=0.0;
1059 value=solution1[i];
1060 if ( fabs(value) > zeroTolerance_ ){
1061 region2[numberNonZero2] = value;
1062 regionIndex2[numberNonZero2++]=i;
1063 }
1064 }
1065 }
1066 regionSparse2->setNumElements(numberNonZero2);
1067 //
1068 numberNonZero3=0;
1069 if (!regionSparse3->packedMode()) {
1070 double value;
1071 for (int i=0;i<numberRows_;i++) {
1072 value=solution2[i];
1073 if ( fabs(value) > zeroTolerance_ ){
1074 vec2[i]=value;
1075 regionIndex3[numberNonZero3++]=i;
1076 }
1077 else
1078 vec2[i]=0.0;
1079 }
1080 }
1081 else { // packed mode
1082 double value;
1083 for (int i=0;i<numberRows_;i++) {
1084 value=solution2[i];
1085 if ( fabs(value) > zeroTolerance_ ){
1086 region3[numberNonZero3] = value;
1087 regionIndex3[numberNonZero3++]=i;
1088 }
1089 }
1090 }
1091 regionSparse3->setNumElements(numberNonZero3);
1092 return 0;
1093 }
1094
1095
1096
1097 int
updateColumnTranspose(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2) const1098 CoinSimpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
1099 CoinIndexedVector * regionSparse2) const
1100 {
1101 upColumnTranspose(regionSparse, regionSparse2);
1102 return 0;
1103 }
1104
1105 int
upColumnTranspose(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2) const1106 CoinSimpFactorization::upColumnTranspose ( CoinIndexedVector * regionSparse,
1107 CoinIndexedVector * regionSparse2) const
1108 {
1109 assert (numberRows_==numberColumns_);
1110 double *region2 = regionSparse2->denseVector ( );
1111 int *regionIndex = regionSparse2->getIndices ( );
1112 int numberNonZero = regionSparse2->getNumElements ( );
1113 double *region = regionSparse->denseVector ( );
1114 if (!regionSparse2->packedMode()) {
1115 region=regionSparse2->denseVector ( );
1116 }
1117 else { // packed
1118 for (int j=0;j<numberNonZero;j++) {
1119 region[regionIndex[j]]=region2[j];
1120 region2[j]=0.0;
1121 }
1122 }
1123 double *solution=workArea2_;
1124 btran(region, solution);
1125 // get nonzeros
1126 numberNonZero=0;
1127 if (!regionSparse2->packedMode()) {
1128 double value;
1129 for (int i=0;i<numberRows_;i++) {
1130 value=solution[i];
1131 if ( fabs(value) > zeroTolerance_ ){
1132 region[i]=value;
1133 regionIndex[numberNonZero++]=i;
1134 }
1135 else
1136 region[i]=0.0;
1137 }
1138 }
1139 else { // packed mode
1140 memset(region,0,numberRows_*sizeof(double));
1141 for (int i=0;i<numberRows_;i++) {
1142 const double value=solution[i];
1143 if ( fabs(value) > zeroTolerance_ ){
1144 region2[numberNonZero] = value;
1145 regionIndex[numberNonZero++]=i;
1146 }
1147 }
1148 }
1149 regionSparse2->setNumElements(numberNonZero);
1150 return 0;
1151 }
1152
1153
1154
1155 int
mainLoopFactor(FactorPointers & pointers)1156 CoinSimpFactorization::mainLoopFactor (FactorPointers &pointers )
1157 {
1158 numberGoodU_=0;
1159 numberSlacks_=0;
1160 bool ifSlack=true;
1161 for ( int i=0; i<numberColumns_; ++i ){
1162 int r, s;
1163 //s=i;
1164 if ( findPivot(pointers,r,s,ifSlack) ){
1165 return -1;
1166 }
1167 if ( ifSlack ) ++numberSlacks_;
1168 const int rowPos=rowPosition_[r];
1169 const int colPos=colPosition_[s];
1170 assert( i <= rowPos && rowPos < numberRows_);
1171 assert( i <= colPos && colPos < numberColumns_);
1172 // permute columns
1173 int j=colOfU_[i];
1174 colOfU_[i]=colOfU_[colPos];
1175 colOfU_[colPos]=j;
1176 colPosition_[colOfU_[i]]=i;
1177 colPosition_[colOfU_[colPos]]=colPos;
1178 // permute rows
1179 j=rowOfU_[i];
1180 rowOfU_[i]=rowOfU_[rowPos];
1181 rowOfU_[rowPos]=j;
1182 rowPosition_[rowOfU_[i]]=i;
1183 rowPosition_[rowOfU_[rowPos]]=rowPos;
1184 GaussEliminate(pointers,r,s);
1185 //if ( maxU_ > maxGrowth_ * maxA_ ){
1186 // return -3;
1187 //}
1188 ++numberGoodU_;
1189 }
1190 return 0;
1191 }
1192 /// find a pivot in the active part of U
findPivot(FactorPointers & pointers,int & r,int & s,bool & ifSlack)1193 int CoinSimpFactorization::findPivot(FactorPointers &pointers, int &r,
1194 int &s, bool &ifSlack){
1195 int *firstRowKnonzeros=pointers.firstRowKnonzeros;
1196 int *nextRow=pointers.nextRow;
1197 int *firstColKnonzeros=pointers.firstColKnonzeros;
1198 int *prevColumn=pointers.prevColumn;
1199 int *nextColumn=pointers.nextColumn;
1200 r=s=-1;
1201 int numCandidates=0;
1202 double bestMarkowitzCount=COIN_DBL_MAX;
1203 // if there is a column with one element choose it as pivot
1204 int column=firstColKnonzeros[1];
1205 if ( column!=-1 ){
1206 assert( UcolLengths_[column] == 1 );
1207 r=UcolInd_[UcolStarts_[column]];
1208 s=column;
1209 if ( !colSlack_[column] )
1210 ifSlack=false;
1211 return 0;
1212 }
1213 // from now on no more slacks
1214 ifSlack=false;
1215 // if there is a row with one element choose it
1216 int row=firstRowKnonzeros[1];
1217 if ( row!=-1 ){
1218 assert( UrowLengths_[row] == 1 );
1219 s=UrowInd_[UrowStarts_[row]];
1220 r=row;
1221 return 0;
1222 }
1223 // consider other rows and columns
1224 for ( int length=2; length <=numberRows_; ++length){
1225 int nextCol=-1;
1226 for ( column=firstColKnonzeros[length]; column!=-1; column=nextCol ){
1227 nextCol=nextColumn[column];
1228 int minRow, minRowLength;
1229 int rc=findShortRow(column, length, minRow, minRowLength, pointers);
1230 if ( rc== 0 ){
1231 r=minRow;
1232 s=column;
1233 return 0;
1234 }
1235 if ( minRow != -1 ){
1236 ++numCandidates;
1237 double MarkowitzCount=static_cast<double>(minRowLength-1)*(length-1);
1238 if ( MarkowitzCount < bestMarkowitzCount ){
1239 r=minRow; s=column;
1240 bestMarkowitzCount=MarkowitzCount;
1241 }
1242 if ( numCandidates == pivotCandLimit_ ) return 0;
1243 }
1244 else {
1245 if ( doSuhlHeuristic_ ){
1246 // this column did not give a candidate, it will be
1247 // removed until it becomes a singleton
1248 removeColumnFromActSet(column, pointers);
1249 prevColumn[column]=nextColumn[column]=column;
1250 }
1251 }
1252 } // end for ( column= ....
1253 // now rows
1254 for ( row=firstRowKnonzeros[length]; row!=-1; row=nextRow[row] ){
1255 int minCol, minColLength;
1256 int rc=findShortColumn(row, length, minCol, minColLength, pointers);
1257 if ( rc==0 ){
1258 r=row;
1259 s=minCol;
1260 return 0;
1261 }
1262 if ( minCol != -1 ){
1263 ++numCandidates;
1264 double MarkowitzCount=static_cast<double>(minColLength-1)*(length-1);
1265 if ( MarkowitzCount < bestMarkowitzCount ){
1266 r=row; s=minCol;
1267 bestMarkowitzCount=MarkowitzCount;
1268 }
1269 if ( numCandidates == pivotCandLimit_ ) return 0;
1270 }
1271 //else abort();
1272 }// end for ( row= ...
1273 }// end for ( int length= ...
1274 if ( r== -1 || s==-1 ) return 1;
1275 else return 0;
1276 }
1277 //
findPivotShCol(FactorPointers & pointers,int & r,int & s)1278 int CoinSimpFactorization::findPivotShCol(FactorPointers &pointers, int &r, int &s)
1279 {
1280 int *firstColKnonzeros=pointers.firstColKnonzeros;
1281 r=s=-1;
1282 // if there is a column with one element choose it as pivot
1283 int column=firstColKnonzeros[1];
1284 if ( column!=-1 ){
1285 assert( UcolLengths_[column] == 1 );
1286 r=UcolInd_[UcolStarts_[column]];
1287 s=column;
1288 return 0;
1289 }
1290 // consider other columns
1291 for ( int length=2; length <=numberRows_; ++length){
1292 column=firstColKnonzeros[length];
1293 if ( column != -1 ) break;
1294 }
1295 if ( column == -1 ) return 1;
1296 // find largest element
1297 const int colBeg=UcolStarts_[column];
1298 const int colEnd=colBeg+UcolLengths_[column];
1299 double largest=0.0;
1300 int rowLargest=-1;
1301 for ( int j=colBeg; j<colEnd; ++j ){
1302 const int row=UcolInd_[j];
1303 const int columnIndx=findInRow(row,column);
1304 assert(columnIndx!=-1);
1305 double coeff=fabs(Urows_[columnIndx]);
1306 if ( coeff < largest ) continue;
1307 largest=coeff;
1308 rowLargest=row;
1309 }
1310 assert(rowLargest != -1);
1311 s=column;
1312 r=rowLargest;
1313 return 0;
1314 }
1315
findPivotSimp(FactorPointers &,int & r,int & s)1316 int CoinSimpFactorization::findPivotSimp(FactorPointers &, int &r, int &s){
1317 r=-1;
1318 int column=s;
1319 const int colBeg=UcolStarts_[column];
1320 const int colEnd=colBeg+UcolLengths_[column];
1321 double largest=0.0;
1322 int rowLargest=-1;
1323 for ( int j=colBeg; j<colEnd; ++j ){
1324 const int row=UcolInd_[j];
1325 const int columnIndx=findInRow(row,column);
1326 assert(columnIndx!=-1);
1327 double coeff=fabs(Urows_[columnIndx]);
1328 if ( coeff < largest ) continue;
1329 largest=coeff;
1330 rowLargest=row;
1331 }
1332 if ( rowLargest != -1 ){
1333 r=rowLargest;
1334 return 0;
1335 }
1336 else return 1;
1337 }
1338
1339
findShortRow(const int column,const int length,int & minRow,int & minRowLength,FactorPointers & pointers)1340 int CoinSimpFactorization::findShortRow(const int column,
1341 const int length,
1342 int &minRow,
1343 int &minRowLength,
1344 FactorPointers &pointers)
1345 {
1346 const int colBeg=UcolStarts_[column];
1347 const int colEnd=colBeg+UcolLengths_[column];
1348 minRow=-1;
1349 minRowLength= COIN_INT_MAX;
1350 for ( int j=colBeg; j<colEnd; ++j ){
1351 int row=UcolInd_[j];
1352 if ( UrowLengths_[row] >= minRowLength ) continue;
1353 double largestInRow=findMaxInRrow(row,pointers);
1354 // find column in row
1355 int columnIndx=findInRow(row,column);
1356 assert(columnIndx!=-1);
1357 double coeff=Urows_[columnIndx];
1358 if ( fabs(coeff) < pivotTolerance_ * largestInRow ) continue;
1359 minRow=row; minRowLength=UrowLengths_[row];
1360 if ( UrowLengths_[row] <= length ) return 0;
1361 }
1362 return 1;
1363 }
findShortColumn(const int row,const int length,int & minCol,int & minColLength,FactorPointers & pointers)1364 int CoinSimpFactorization::findShortColumn(const int row,
1365 const int length,
1366 int &minCol,
1367 int &minColLength,
1368 FactorPointers &pointers)
1369 {
1370 const int rowBeg=UrowStarts_[row];
1371 const int rowEnd=rowBeg+UrowLengths_[row];
1372 minCol=-1;
1373 minColLength=COIN_INT_MAX;
1374 double largestInRow=findMaxInRrow(row,pointers);
1375 for ( int i=rowBeg; i<rowEnd; ++i ){
1376 int column=UrowInd_[i];
1377 if ( UcolLengths_[column] >= minColLength ) continue;
1378 double coeff=Urows_[i];
1379 if ( fabs(coeff) < pivotTolerance_ * largestInRow ) continue;
1380 minCol=column;
1381 minColLength=UcolLengths_[column];
1382 if ( minColLength <= length ) return 0;
1383 }
1384 return 1;
1385 }
1386 // Gaussian elimination
GaussEliminate(FactorPointers & pointers,int & pivotRow,int & pivotCol)1387 void CoinSimpFactorization::GaussEliminate(FactorPointers &pointers, int &pivotRow, int &pivotCol)
1388 {
1389 assert( pivotRow >= 0 && pivotRow < numberRows_ );
1390 assert( pivotCol >= 0 && pivotCol < numberRows_ );
1391 int *firstColKnonzeros=pointers.firstColKnonzeros;
1392 int *prevColumn=pointers.prevColumn;
1393 int *nextColumn=pointers.nextColumn;
1394 int *colLabels=vecLabels_;
1395 double *denseRow=denseVector_;
1396 removeRowFromActSet(pivotRow, pointers);
1397 removeColumnFromActSet(pivotCol, pointers);
1398 // find column s
1399 int indxColS=findInRow(pivotRow, pivotCol);
1400 assert( indxColS >= 0 );
1401 // store the inverse of the pivot and remove it from row
1402 double invPivot=1.0/Urows_[indxColS];
1403 invOfPivots_[pivotRow]=invPivot;
1404 int rowBeg=UrowStarts_[pivotRow];
1405 int rowEnd=rowBeg+UrowLengths_[pivotRow];
1406 Urows_[indxColS]=Urows_[rowEnd-1];
1407 UrowInd_[indxColS]=UrowInd_[rowEnd-1];
1408 --UrowLengths_[pivotRow];
1409 --rowEnd;
1410 // now remove pivot from column
1411 int indxRowR=findInColumn(pivotCol,pivotRow);
1412 assert( indxRowR >= 0 );
1413 const int pivColEnd=UcolStarts_[pivotCol]+UcolLengths_[pivotCol];
1414 UcolInd_[indxRowR]=UcolInd_[pivColEnd-1];
1415 --UcolLengths_[pivotCol];
1416 // go through pivot row
1417 for ( int i=rowBeg; i<rowEnd; ++i ){
1418 int column=UrowInd_[i];
1419 colLabels[column]=1;
1420 denseRow[column]=Urows_[i];
1421 // remove this column from bucket because it will change
1422 removeColumnFromActSet(column, pointers);
1423 // remove element (pivotRow, column) from column
1424 int indxRow=findInColumn(column,pivotRow);
1425 assert( indxRow>=0 );
1426 const int colEnd=UcolStarts_[column]+UcolLengths_[column];
1427 UcolInd_[indxRow]=UcolInd_[colEnd-1];
1428 --UcolLengths_[column];
1429 }
1430 //
1431 pivoting(pivotRow, pivotCol, invPivot, pointers);
1432 //
1433 rowBeg=UrowStarts_[pivotRow];
1434 rowEnd=rowBeg+UrowLengths_[pivotRow];
1435 for ( int i=rowBeg; i<rowEnd; ++i ){
1436 int column=UrowInd_[i];
1437 // clean back these two arrays
1438 colLabels[column]=0;
1439 denseRow[column]=0.0;
1440 // column goes into a bucket, if Suhl' heuristic had removed it, it
1441 // can go back only if it is a singleton
1442 if ( UcolLengths_[column] != 1 ||
1443 prevColumn[column]!=column || nextColumn[column]!=column )
1444 {
1445 prevColumn[column]=-1;
1446 nextColumn[column]=firstColKnonzeros[UcolLengths_[column]];
1447 if ( nextColumn[column] != -1 )
1448 prevColumn[nextColumn[column]]=column;
1449 firstColKnonzeros[UcolLengths_[column]]=column;
1450 }
1451 }
1452 }
pivoting(const int pivotRow,const int pivotColumn,const double invPivot,FactorPointers & pointers)1453 void CoinSimpFactorization::pivoting(const int pivotRow,
1454 const int pivotColumn,
1455 const double invPivot,
1456 FactorPointers &pointers)
1457 {
1458 // initialize the new column of L
1459 LcolStarts_[pivotRow]=LcolSize_;
1460 // go trough pivot column
1461 const int colBeg=UcolStarts_[pivotColumn];
1462 int colEnd=colBeg+UcolLengths_[pivotColumn];
1463 for ( int i=colBeg; i<colEnd; ++i ){
1464 int row=UcolInd_[i];
1465 // remove row from bucket because it will change
1466 removeRowFromActSet(row, pointers);
1467 // find pivot column
1468 int pivotColInRow=findInRow(row,pivotColumn);
1469 assert(pivotColInRow >= 0);
1470 const double multiplier=Urows_[pivotColInRow]*invPivot;
1471 // remove element (row,pivotColumn) from row
1472 const int currentRowEnd=UrowStarts_[row]+UrowLengths_[row];
1473 Urows_[pivotColInRow]=Urows_[currentRowEnd-1];
1474 UrowInd_[pivotColInRow]=UrowInd_[currentRowEnd-1];
1475 --UrowLengths_[row];
1476 int newNonZeros=UrowLengths_[pivotRow];
1477 updateCurrentRow(pivotRow, row, multiplier, pointers,
1478 newNonZeros);
1479 // store multiplier
1480 if ( LcolSize_ == LcolCap_ ) increaseLsize();
1481 Lcolumns_[LcolSize_]=multiplier;
1482 LcolInd_[LcolSize_++]=row;
1483 ++LcolLengths_[pivotRow];
1484 }
1485 // remove elements of pivot column
1486 UcolLengths_[pivotColumn]=0;
1487 // remove pivot column from Ucol_
1488 if ( prevColInU_[pivotColumn]==-1 )
1489 firstColInU_=nextColInU_[pivotColumn];
1490 else{
1491 nextColInU_[prevColInU_[pivotColumn]]=nextColInU_[pivotColumn];
1492 #ifdef COIN_SIMP_CAPACITY
1493 UcolCapacities_[prevColInU_[pivotColumn]]+=UcolCapacities_[pivotColumn];
1494 UcolCapacities_[pivotColumn]=0;
1495 #endif
1496 }
1497 if ( nextColInU_[pivotColumn] == -1 )
1498 lastColInU_=prevColInU_[pivotColumn];
1499 else
1500 prevColInU_[nextColInU_[pivotColumn]]=prevColInU_[pivotColumn];
1501 }
1502
1503
1504
updateCurrentRow(const int pivotRow,const int row,const double multiplier,FactorPointers & pointers,int & newNonZeros)1505 void CoinSimpFactorization::updateCurrentRow(const int pivotRow,
1506 const int row,
1507 const double multiplier,
1508 FactorPointers &pointers,
1509 int &newNonZeros)
1510 {
1511 double *rowMax=pointers.rowMax;
1512 int *firstRowKnonzeros=pointers.firstRowKnonzeros;
1513 int *prevRow=pointers.prevRow;
1514 int *nextRow=pointers.nextRow;
1515 int *colLabels=vecLabels_;
1516 double *denseRow=denseVector_;
1517 const int rowBeg=UrowStarts_[row];
1518 int rowEnd=rowBeg+UrowLengths_[row];
1519 // treat old nonzeros
1520 for ( int i=rowBeg; i<rowEnd; ++i ){
1521 const int column=UrowInd_[i];
1522 if ( colLabels[column] ){
1523 Urows_[i]-= multiplier*denseRow[column];
1524 const double absNewCoeff=fabs(Urows_[i]);
1525 colLabels[column]=0;
1526 --newNonZeros;
1527 if ( absNewCoeff < zeroTolerance_ ){
1528 // remove it from row
1529 UrowInd_[i]=UrowInd_[rowEnd-1];
1530 Urows_[i]=Urows_[rowEnd-1];
1531 --UrowLengths_[row];
1532 --i;
1533 --rowEnd;
1534 // remove it from column
1535 int indxRow=findInColumn(column, row);
1536 assert( indxRow >= 0 );
1537 const int colEnd=UcolStarts_[column]+UcolLengths_[column];
1538 UcolInd_[indxRow]=UcolInd_[colEnd-1];
1539 --UcolLengths_[column];
1540 }
1541 else
1542 {
1543 if ( maxU_ < absNewCoeff )
1544 maxU_=absNewCoeff;
1545 }
1546 }
1547 }
1548 // now add the new nonzeros to the row
1549 #ifdef COIN_SIMP_CAPACITY
1550 if ( UrowLengths_[row] + newNonZeros > UrowCapacities_[row] )
1551 increaseRowSize(row, UrowLengths_[row] + newNonZeros);
1552 #endif
1553 const int pivotRowBeg=UrowStarts_[pivotRow];
1554 const int pivotRowEnd=pivotRowBeg+UrowLengths_[pivotRow];
1555 int numNew=0;
1556 int *newCols=pointers.newCols;
1557 for ( int i=pivotRowBeg; i<pivotRowEnd; ++i ){
1558 const int column=UrowInd_[i];
1559 if ( colLabels[column] ){
1560 const double value= -multiplier*denseRow[column];
1561 const double absNewCoeff=fabs(value);
1562 if ( absNewCoeff >= zeroTolerance_ ){
1563 const int newInd=UrowStarts_[row]+UrowLengths_[row];
1564 Urows_[newInd]=value;
1565 UrowInd_[newInd]=column;
1566 ++UrowLengths_[row];
1567 newCols[numNew++]=column;
1568 if ( maxU_ < absNewCoeff ) maxU_=absNewCoeff;
1569 }
1570 }
1571 else colLabels[column]=1;
1572 }
1573 // add the new nonzeros to the columns
1574 for ( int i=0; i<numNew; ++i){
1575 const int column=newCols[i];
1576 #ifdef COIN_SIMP_CAPACITY
1577 if ( UcolLengths_[column] + 1 > UcolCapacities_[column] ){
1578 increaseColSize(column, UcolLengths_[column] + 1, false);
1579 }
1580 #endif
1581 const int newInd=UcolStarts_[column]+UcolLengths_[column];
1582 UcolInd_[newInd]=row;
1583 ++UcolLengths_[column];
1584 }
1585 // the row goes to a new bucket
1586 prevRow[row]=-1;
1587 nextRow[row]=firstRowKnonzeros[UrowLengths_[row]];
1588 if ( nextRow[row]!=-1 ) prevRow[nextRow[row]]=row;
1589 firstRowKnonzeros[UrowLengths_[row]]=row;
1590 //
1591 rowMax[row]=-1.0;
1592 }
1593 #ifdef COIN_SIMP_CAPACITY
1594
increaseRowSize(const int row,const int newSize)1595 void CoinSimpFactorization::increaseRowSize(const int row,
1596 const int newSize)
1597 {
1598 assert( newSize > UrowCapacities_[row] );
1599 const int newNumElements=newSize + minIncrease_;
1600 if ( UrowMaxCap_ < UrowEnd_ + newNumElements ){
1601 enlargeUrow( UrowEnd_ + newNumElements - UrowMaxCap_ );
1602 }
1603 int currentCapacity=UrowCapacities_[row];
1604 memcpy(&Urows_[UrowEnd_],&Urows_[UrowStarts_[row]],
1605 UrowLengths_[row] * sizeof(double));
1606 memcpy(&UrowInd_[UrowEnd_],&UrowInd_[UrowStarts_[row]],
1607 UrowLengths_[row] * sizeof(int));
1608 UrowStarts_[row]=UrowEnd_;
1609 UrowCapacities_[row]=newNumElements;
1610 UrowEnd_+=UrowCapacities_[row];
1611 if ( firstRowInU_==lastRowInU_ ) return; // only one element
1612 // remove row from list
1613 if( prevRowInU_[row]== -1)
1614 firstRowInU_=nextRowInU_[row];
1615 else {
1616 nextRowInU_[prevRowInU_[row]]=nextRowInU_[row];
1617 UrowCapacities_[prevRowInU_[row]]+=currentCapacity;
1618 }
1619 if ( nextRowInU_[row]==-1 )
1620 lastRowInU_=prevRowInU_[row];
1621 else
1622 prevRowInU_[nextRowInU_[row]]=prevRowInU_[row];
1623 // add row at the end of list
1624 nextRowInU_[lastRowInU_]=row;
1625 nextRowInU_[row]=-1;
1626 prevRowInU_[row]=lastRowInU_;
1627 lastRowInU_=row;
1628 }
1629 #endif
1630 #ifdef COIN_SIMP_CAPACITY
increaseColSize(const int column,const int newSize,const bool ifElements)1631 void CoinSimpFactorization::increaseColSize(const int column,
1632 const int newSize,
1633 const bool ifElements)
1634 {
1635 assert( newSize > UcolCapacities_[column] );
1636 const int newNumElements=newSize+minIncrease_;
1637 if ( UcolMaxCap_ < UcolEnd_ + newNumElements ){
1638 enlargeUcol(UcolEnd_ + newNumElements - UcolMaxCap_,
1639 ifElements);
1640 }
1641 int currentCapacity=UcolCapacities_[column];
1642 memcpy(&UcolInd_[UcolEnd_], &UcolInd_[UcolStarts_[column]],
1643 UcolLengths_[column] * sizeof(int));
1644 if ( ifElements ){
1645 memcpy(&Ucolumns_[UcolEnd_], &Ucolumns_[UcolStarts_[column]],
1646 UcolLengths_[column] * sizeof(double) );
1647 }
1648 UcolStarts_[column]=UcolEnd_;
1649 UcolCapacities_[column]=newNumElements;
1650 UcolEnd_+=UcolCapacities_[column];
1651 if ( firstColInU_==lastColInU_ ) return; // only one column
1652 // remove from list
1653 if ( prevColInU_[column]==-1 )
1654 firstColInU_=nextColInU_[column];
1655 else {
1656 nextColInU_[prevColInU_[column]]=nextColInU_[column];
1657 UcolCapacities_[prevColInU_[column]]+=currentCapacity;
1658 }
1659 if ( nextColInU_[column]==-1 )
1660 lastColInU_=prevColInU_[column];
1661 else
1662 prevColInU_[nextColInU_[column]]=prevColInU_[column];
1663 // add column at the end
1664 nextColInU_[lastColInU_]=column;
1665 nextColInU_[column]=-1;
1666 prevColInU_[column]=lastColInU_;
1667 lastColInU_=column;
1668 }
1669 #endif
findMaxInRrow(const int row,FactorPointers & pointers)1670 double CoinSimpFactorization::findMaxInRrow(const int row,
1671 FactorPointers &pointers)
1672 {
1673 double *rowMax=pointers.rowMax;
1674 double largest=rowMax[row];
1675 if ( largest >= 0.0 ) return largest;
1676 const int rowBeg=UrowStarts_[row];
1677 const int rowEnd=rowBeg+UrowLengths_[row];
1678 for ( int i=rowBeg; i<rowEnd; ++i ) {
1679 const double absValue=fabs(Urows_[i]);
1680 if ( absValue > largest )
1681 largest=absValue;
1682 }
1683 rowMax[row]=largest;
1684 return largest;
1685 }
increaseLsize()1686 void CoinSimpFactorization::increaseLsize()
1687 {
1688 int newcap= LcolCap_ + minIncrease_;
1689
1690 double *aux=new double[newcap];
1691 memcpy(aux, Lcolumns_, LcolCap_ * sizeof(double));
1692 delete [] Lcolumns_;
1693 Lcolumns_=aux;
1694
1695 int *iaux=new int[newcap];
1696 memcpy(iaux, LcolInd_, LcolCap_ * sizeof(int));
1697 delete [] LcolInd_;
1698 LcolInd_=iaux;
1699
1700 LcolCap_=newcap;
1701 }
enlargeUcol(const int numNewElements,const bool ifElements)1702 void CoinSimpFactorization::enlargeUcol(const int numNewElements, const bool ifElements)
1703 {
1704 int *iaux=new int[UcolMaxCap_+numNewElements];
1705 memcpy(iaux, UcolInd_, UcolMaxCap_*sizeof(int) );
1706 delete [] UcolInd_;
1707 UcolInd_=iaux;
1708
1709 if ( ifElements ){
1710 double *aux=new double[UcolMaxCap_+numNewElements];
1711 memcpy(aux, Ucolumns_, UcolMaxCap_*sizeof(double) );
1712 delete [] Ucolumns_;
1713 Ucolumns_=aux;
1714 }
1715
1716 UcolMaxCap_+=numNewElements;
1717 }
enlargeUrow(const int numNewElements)1718 void CoinSimpFactorization::enlargeUrow(const int numNewElements)
1719 {
1720 int *iaux=new int[UrowMaxCap_+numNewElements];
1721 memcpy(iaux, UrowInd_, UrowMaxCap_*sizeof(int) );
1722 delete [] UrowInd_;
1723 UrowInd_=iaux;
1724
1725 double *aux=new double[UrowMaxCap_+numNewElements];
1726 memcpy(aux, Urows_, UrowMaxCap_*sizeof(double) );
1727 delete [] Urows_;
1728 Urows_=aux;
1729
1730 UrowMaxCap_+=numNewElements;
1731 }
1732
copyUbyColumns()1733 void CoinSimpFactorization::copyUbyColumns()
1734 {
1735 memset(UcolLengths_,0,numberColumns_*sizeof(int));
1736 for ( int column=0; column<numberColumns_; ++column ){
1737 prevColInU_[column]=column-1;
1738 nextColInU_[column]=column+1;
1739 }
1740 nextColInU_[numberColumns_-1]=-1;
1741 firstColInU_=0;
1742 lastColInU_=numberColumns_-1;
1743 //int nonZeros=0;
1744 //for ( int row=0; row<numberRows_; ++row ){
1745 //const int rowBeg=UrowStarts_[row];
1746 //const int rowEnd=rowBeg+UrowLengths_[row];
1747 //for ( int j=rowBeg; j<rowEnd; ++j )
1748 // ++UcolCapacities_[UrowInd_[j]];
1749 // nonZeros+=UrowLengths_[row];
1750 // }
1751 //
1752 //memset(UcolCapacities_, numberRows_, numberColumns_ * sizeof(int));
1753 #ifdef COIN_SIMP_CAPACITY
1754 for ( int i=0; i<numberColumns_; ++i ) UcolCapacities_[i]=numberRows_;
1755 #endif
1756 int k=0;
1757 for ( int column=0; column<numberColumns_; ++column ){
1758 UcolStarts_[column]=k;
1759 //UcolCapacities_[column]+=minIncrease_;
1760 //k+=UcolCapacities_[column];
1761 k+=numberRows_;
1762 }
1763 UcolEnd_=k;
1764 // go through the rows and fill the columns, assume UcolLengths_[]=0
1765 for ( int row=0; row<numberRows_; ++row ){
1766 const int rowBeg=UrowStarts_[row];
1767 int rowEnd=rowBeg+UrowLengths_[row];
1768 for ( int j=rowBeg; j<rowEnd; ++j ){
1769 // remove els close to zero
1770 while( fabs( Urows_[j] ) < zeroTolerance_ ){
1771 --UrowLengths_[row];
1772 --rowEnd;
1773 if ( j < rowEnd ){
1774 Urows_[j]=Urows_[rowEnd];
1775 UrowInd_[j]=UrowInd_[rowEnd];
1776 }
1777 else break;
1778 }
1779 if ( j==rowEnd ) continue;
1780 //
1781 const int column=UrowInd_[j];
1782 const int indx=UcolStarts_[column]+UcolLengths_[column];
1783 Ucolumns_[indx]=Urows_[j];
1784 UcolInd_[indx]=row;
1785 ++UcolLengths_[column];
1786 }
1787 }
1788 }
copyLbyRows()1789 void CoinSimpFactorization::copyLbyRows()
1790 {
1791 int nonZeros=0;
1792 memset(LrowLengths_,0,numberRows_*sizeof(int));
1793 for ( int column=0; column<numberRows_; ++column ){
1794 const int colBeg=LcolStarts_[column];
1795 const int colEnd=colBeg+LcolLengths_[column];
1796 for ( int j=colBeg; j<colEnd; ++j )
1797 ++LrowLengths_[LcolInd_[j]];
1798 nonZeros+=LcolLengths_[column];
1799 }
1800 //
1801 LrowSize_=nonZeros;
1802 int k=0;
1803 for ( int row=0; row<numberRows_; ++row ){
1804 LrowStarts_[row]=k;
1805 k+=LrowLengths_[row];
1806 }
1807 #ifdef COIN_SIMP_CAPACITY
1808 //memset(LrowLengths_,0,numberRows_*sizeof(int));
1809 // fill the rows
1810 for ( int column=0; column<numberRows_; ++column ){
1811 const int colBeg=LcolStarts_[column];
1812 const int colEnd=colBeg+LcolLengths_[column];
1813 for ( int j=colBeg; j<colEnd; ++j ){
1814 const int row=LcolInd_[j];
1815 const int indx=LrowStarts_[row]++;
1816 Lrows_[indx]=Lcolumns_[j];
1817 LrowInd_[indx]=column;
1818 }
1819 }
1820 // Put back starts
1821 k=0;
1822 for ( int row=0; row<numberRows_; ++row ){
1823 int next = LrowStarts_[row];
1824 LrowStarts_[row]=k;
1825 k=next;
1826 }
1827 #else
1828 memset(LrowLengths_,0,numberRows_*sizeof(int));
1829 // fill the rows
1830 for ( int column=0; column<numberRows_; ++column ){
1831 const int colBeg=LcolStarts_[column];
1832 const int colEnd=colBeg+LcolLengths_[column];
1833 for ( int j=colBeg; j<colEnd; ++j ){
1834 const int row=LcolInd_[j];
1835 const int indx=LrowStarts_[row] + LrowLengths_[row];
1836 Lrows_[indx]=Lcolumns_[j];
1837 LrowInd_[indx]=column;
1838 ++LrowLengths_[row];
1839 }
1840 }
1841 #endif
1842 }
1843
1844
findInRow(const int row,const int column)1845 int CoinSimpFactorization::findInRow(const int row,
1846 const int column)
1847 {
1848 const int rowBeg=UrowStarts_[row];
1849 const int rowEnd=rowBeg+UrowLengths_[row];
1850 int columnIndx=-1;
1851 for ( int i=rowBeg; i<rowEnd; ++i ){
1852 if ( UrowInd_[i]==column ){
1853 columnIndx=i;
1854 break;
1855 }
1856 }
1857 return columnIndx;
1858 }
findInColumn(const int column,const int row)1859 int CoinSimpFactorization::findInColumn(const int column,
1860 const int row)
1861 {
1862 int indxRow=-1;
1863 const int colBeg=UcolStarts_[column];
1864 const int colEnd=colBeg+UcolLengths_[column];
1865 for ( int i=colBeg; i<colEnd; ++i ){
1866 if (UcolInd_[i]==row){
1867 indxRow=i;
1868 break;
1869 }
1870 }
1871 return indxRow;
1872 }
removeRowFromActSet(const int row,FactorPointers & pointers)1873 void CoinSimpFactorization::removeRowFromActSet(const int row,
1874 FactorPointers &pointers)
1875 {
1876 int *firstRowKnonzeros=pointers.firstRowKnonzeros;
1877 int *prevRow=pointers.prevRow;
1878 int *nextRow=pointers.nextRow;
1879 if ( prevRow[row]==-1 )
1880 firstRowKnonzeros[UrowLengths_[row]]=nextRow[row];
1881 else nextRow[prevRow[row]]=nextRow[row];
1882 if ( nextRow[row] != -1)
1883 prevRow[nextRow[row]]=prevRow[row];
1884 }
removeColumnFromActSet(const int column,FactorPointers & pointers)1885 void CoinSimpFactorization::removeColumnFromActSet(const int column,
1886 FactorPointers &pointers)
1887 {
1888 int *firstColKnonzeros=pointers.firstColKnonzeros;
1889 int *prevColumn=pointers.prevColumn;
1890 int *nextColumn=pointers.nextColumn;
1891 if ( prevColumn[column]==-1 )
1892 firstColKnonzeros[UcolLengths_[column]]=nextColumn[column];
1893 else nextColumn[prevColumn[column]]=nextColumn[column];
1894 if ( nextColumn[column] != -1 )
1895 prevColumn[nextColumn[column]]=prevColumn[column];
1896 }
allocateSomeArrays()1897 void CoinSimpFactorization::allocateSomeArrays()
1898 {
1899 if (denseVector_) delete [] denseVector_;
1900 denseVector_ = new double[numberRows_];
1901 memset(denseVector_,0,numberRows_*sizeof(double));
1902 if(workArea2_) delete [] workArea2_;
1903 workArea2_ = new double[numberRows_];
1904 if(workArea3_) delete [] workArea3_;
1905 workArea3_ = new double[numberRows_];
1906
1907 if(vecLabels_) delete [] vecLabels_;
1908 vecLabels_ = new int[numberRows_];
1909 memset(vecLabels_,0, numberRows_*sizeof(int));
1910 if(indVector_) delete [] indVector_;
1911 indVector_ = new int[numberRows_];
1912
1913 if (auxVector_) delete [] auxVector_;
1914 auxVector_ = new double[numberRows_];
1915 if (auxInd_) delete [] auxInd_;
1916 auxInd_ = new int[numberRows_];
1917
1918 if (vecKeep_) delete [] vecKeep_;
1919 vecKeep_ = new double[numberRows_];
1920 if (indKeep_) delete [] indKeep_;
1921 indKeep_ = new int[numberRows_];
1922
1923 if (LrowStarts_) delete [] LrowStarts_;
1924 LrowStarts_ = new int[numberRows_];
1925 if (LrowLengths_) delete [] LrowLengths_;
1926 LrowLengths_ = new int[numberRows_];
1927
1928 LrowCap_=(numberRows_*(numberRows_-1))/2;
1929 if (Lrows_) delete [] Lrows_;
1930 Lrows_ = new double[LrowCap_];
1931 if (LrowInd_) delete [] LrowInd_;
1932 LrowInd_ = new int[LrowCap_];
1933
1934 if (LcolStarts_) delete [] LcolStarts_;
1935 LcolStarts_ = new int[numberRows_];
1936 if (LcolLengths_) delete [] LcolLengths_;
1937 LcolLengths_ = new int[numberRows_];
1938 LcolCap_=LrowCap_;
1939 if (Lcolumns_) delete [] Lcolumns_;
1940 Lcolumns_ = new double[LcolCap_];
1941 if (LcolInd_) delete [] LcolInd_;
1942 LcolInd_ = new int[LcolCap_];
1943
1944 if (UrowStarts_) delete [] UrowStarts_;
1945 UrowStarts_ = new int[numberRows_];
1946 if (UrowLengths_) delete [] UrowLengths_;
1947 UrowLengths_ = new int[numberRows_];
1948 #ifdef COIN_SIMP_CAPACITY
1949 if (UrowCapacities_) delete [] UrowCapacities_;
1950 UrowCapacities_ = new int[numberRows_];
1951 #endif
1952
1953 minIncrease_=10;
1954 UrowMaxCap_=numberRows_*(numberRows_+minIncrease_);
1955 if (Urows_) delete [] Urows_;
1956 Urows_ = new double[UrowMaxCap_];
1957 if (UrowInd_) delete [] UrowInd_;
1958 UrowInd_ = new int[UrowMaxCap_];
1959
1960 if (prevRowInU_) delete [] prevRowInU_;
1961 prevRowInU_ = new int[numberRows_];
1962 if (nextRowInU_) delete [] nextRowInU_;
1963 nextRowInU_ = new int[numberRows_];
1964
1965 if (UcolStarts_) delete [] UcolStarts_;
1966 UcolStarts_ = new int[numberRows_];
1967 if (UcolLengths_) delete [] UcolLengths_;
1968 UcolLengths_ = new int[numberRows_];
1969 #ifdef COIN_SIMP_CAPACITY
1970 if (UcolCapacities_) delete [] UcolCapacities_;
1971 UcolCapacities_ = new int[numberRows_];
1972 #endif
1973
1974 UcolMaxCap_=UrowMaxCap_;
1975 if (Ucolumns_) delete [] Ucolumns_;
1976 Ucolumns_ = new double[UcolMaxCap_];
1977 if (UcolInd_) delete [] UcolInd_;
1978 UcolInd_ = new int[UcolMaxCap_];
1979
1980 if (prevColInU_) delete [] prevColInU_;
1981 prevColInU_ = new int[numberRows_];
1982 if (nextColInU_) delete [] nextColInU_;
1983 nextColInU_ = new int[numberRows_];
1984 if (colSlack_) delete [] colSlack_;
1985 colSlack_ = new int[numberRows_];
1986
1987 if (invOfPivots_) delete [] invOfPivots_;
1988 invOfPivots_ = new double[numberRows_];
1989
1990 if (colOfU_) delete [] colOfU_;
1991 colOfU_ = new int[numberRows_];
1992 if (colPosition_) delete [] colPosition_;
1993 colPosition_ = new int[numberRows_];
1994 if (rowOfU_) delete [] rowOfU_;
1995 rowOfU_ = new int[numberRows_];
1996 if (rowPosition_) delete [] rowPosition_;
1997 rowPosition_ = new int[numberRows_];
1998 if (secRowOfU_) delete [] secRowOfU_;
1999 secRowOfU_ = new int[numberRows_];
2000 if (secRowPosition_) delete [] secRowPosition_;
2001 secRowPosition_ = new int[numberRows_];
2002
2003 if (EtaPosition_) delete [] EtaPosition_;
2004 EtaPosition_ = new int[maximumPivots_];
2005 if (EtaStarts_) delete [] EtaStarts_;
2006 EtaStarts_ = new int[maximumPivots_];
2007 if (EtaLengths_) delete [] EtaLengths_;
2008 EtaLengths_ = new int[maximumPivots_];
2009 maxEtaRows_=maximumPivots_;
2010
2011 EtaMaxCap_=maximumPivots_*minIncrease_;
2012 if (EtaInd_) delete [] EtaInd_;
2013 EtaInd_ = new int[EtaMaxCap_];
2014 if (Eta_) delete [] Eta_;
2015 Eta_ = new double[EtaMaxCap_];
2016
2017 }
2018
Lxeqb(double * b) const2019 void CoinSimpFactorization::Lxeqb(double *b) const
2020 {
2021 double *rhs=b;
2022 int k, colBeg, *ind, *indEnd;
2023 double xk, *Lcol;
2024 // now solve
2025 for ( int j=firstNumberSlacks_; j<numberRows_; ++j ){
2026 k=rowOfU_[j];
2027 xk=rhs[k];
2028 if ( xk!=0.0 ) {
2029 //if ( fabs(xk)>zeroTolerance_ ) {
2030 colBeg=LcolStarts_[k];
2031 ind=LcolInd_+colBeg;
2032 indEnd=ind+LcolLengths_[k];
2033 Lcol=Lcolumns_+colBeg;
2034 for ( ; ind!=indEnd; ++ind ){
2035 rhs[ *ind ]-= (*Lcol) * xk;
2036 ++Lcol;
2037 }
2038 }
2039 }
2040 }
2041
2042
2043
Lxeqb2(double * b1,double * b2) const2044 void CoinSimpFactorization::Lxeqb2(double *b1, double *b2) const
2045 {
2046 double *rhs1=b1;
2047 double *rhs2=b2;
2048 double x1, x2, *Lcol;
2049 int k, colBeg, *ind, *indEnd, j;
2050 // now solve
2051 for ( j=firstNumberSlacks_; j<numberRows_; ++j ){
2052 k=rowOfU_[j];
2053 x1=rhs1[k];
2054 x2=rhs2[k];
2055 if ( x1 == 0.0 ) {
2056 if (x2 == 0.0 ) {
2057 } else {
2058 colBeg=LcolStarts_[k];
2059 ind=LcolInd_+colBeg;
2060 indEnd=ind+LcolLengths_[k];
2061 Lcol=Lcolumns_+colBeg;
2062 for ( ; ind!=indEnd; ++ind ){
2063 #if 0
2064 rhs2[ *ind ]-= (*Lcol) * x2;
2065 #else
2066 double value=rhs2[ *ind ];
2067 rhs2[ *ind ]= value -(*Lcol) * x2;
2068 #endif
2069 ++Lcol;
2070 }
2071 }
2072 } else {
2073 if ( x2 == 0.0 ) {
2074 colBeg=LcolStarts_[k];
2075 ind=LcolInd_+colBeg;
2076 indEnd=ind+LcolLengths_[k];
2077 Lcol=Lcolumns_+colBeg;
2078 for ( ; ind!=indEnd; ++ind ){
2079 #if 0
2080 rhs1[ *ind ]-= (*Lcol) * x1;
2081 #else
2082 double value=rhs1[ *ind ];
2083 rhs1[ *ind ]= value -(*Lcol) * x1;
2084 #endif
2085 ++Lcol;
2086 }
2087 } else {
2088 colBeg=LcolStarts_[k];
2089 ind=LcolInd_+colBeg;
2090 indEnd=ind+LcolLengths_[k];
2091 Lcol=Lcolumns_+colBeg;
2092 for ( ; ind!=indEnd; ++ind ){
2093 #if 0
2094 rhs1[ *ind ]-= (*Lcol) * x1;
2095 rhs2[ *ind ]-= (*Lcol) * x2;
2096 #else
2097 double value1=rhs1[ *ind ];
2098 rhs1[ *ind ]= value1 - (*Lcol) * x1;
2099 double value2=rhs2[ *ind ];
2100 rhs2[ *ind ]= value2 - (*Lcol) * x2;
2101 #endif
2102 ++Lcol;
2103 }
2104 }
2105 }
2106 }
2107 }
2108
Uxeqb(double * b,double * sol) const2109 void CoinSimpFactorization::Uxeqb(double *b, double *sol) const
2110 {
2111 double *rhs=b;
2112 int row, column, colBeg, *ind, *indEnd, k;
2113 double x, *uCol;
2114 // now solve
2115 for ( k=numberRows_-1; k>=numberSlacks_; --k ){
2116 row=secRowOfU_[k];
2117 x=rhs[row];
2118 column=colOfU_[k];
2119 if ( x!=0.0 ) {
2120 //if ( fabs(x) > zeroTolerance_ ) {
2121 x*=invOfPivots_[row];
2122 colBeg=UcolStarts_[column];
2123 ind=UcolInd_+colBeg;
2124 indEnd=ind+UcolLengths_[column];
2125 uCol=Ucolumns_+colBeg;
2126 for ( ; ind!=indEnd; ++ind ){
2127 #if 0
2128 rhs[ *ind ]-= (*uCol) * x;
2129 #else
2130 double value=rhs[ *ind ];
2131 rhs[ *ind ] = value - (*uCol) * x;
2132 #endif
2133 ++uCol;
2134 }
2135 sol[column]=x;
2136 }
2137 else sol[column]=0.0;
2138 }
2139 for ( k=numberSlacks_-1; k>=0; --k ){
2140 row=secRowOfU_[k];
2141 column=colOfU_[k];
2142 sol[column]=-rhs[row];
2143 }
2144 }
2145
2146
2147
Uxeqb2(double * b1,double * sol1,double * b2,double * sol2) const2148 void CoinSimpFactorization::Uxeqb2(double *b1, double *sol1, double *b2, double *sol2) const
2149 {
2150 double *rhs1=b1;
2151 double *rhs2=b2;
2152 int row, column, colBeg, *ind, *indEnd;
2153 double x1, x2, *uCol;
2154 // now solve
2155 for ( int k=numberRows_-1; k>=numberSlacks_; --k ){
2156 row=secRowOfU_[k];
2157 x1=rhs1[row];
2158 x2=rhs2[row];
2159 column=colOfU_[k];
2160 if (x1 == 0.0) {
2161 if (x2 == 0.0) {
2162 sol1[column]=0.0;
2163 sol2[column]=0.0;
2164 } else {
2165 x2*=invOfPivots_[row];
2166 colBeg=UcolStarts_[column];
2167 ind=UcolInd_+colBeg;
2168 indEnd=ind+UcolLengths_[column];
2169 uCol=Ucolumns_+colBeg;
2170 for ( ; ind!=indEnd; ++ind ){
2171 #if 0
2172 rhs2[ *ind ]-= (*uCol) * x2;
2173 #else
2174 double value=rhs2[ *ind ];
2175 rhs2[ *ind ]= value - (*uCol) * x2;
2176 #endif
2177 ++uCol;
2178 }
2179 sol1[column]=0.0;
2180 sol2[column]=x2;
2181 }
2182 } else {
2183 if (x2 == 0.0) {
2184 x1*=invOfPivots_[row];
2185 colBeg=UcolStarts_[column];
2186 ind=UcolInd_+colBeg;
2187 indEnd=ind+UcolLengths_[column];
2188 uCol=Ucolumns_+colBeg;
2189 for ( ; ind!=indEnd; ++ind ){
2190 #if 0
2191 rhs1[ *ind ]-= (*uCol) * x1;
2192 #else
2193 double value=rhs1[ *ind ];
2194 rhs1[ *ind ] = value - (*uCol) * x1;
2195 #endif
2196 ++uCol;
2197 }
2198 sol1[column]=x1;
2199 sol2[column]=0.0;
2200 } else {
2201 x1*=invOfPivots_[row];
2202 x2*=invOfPivots_[row];
2203 colBeg=UcolStarts_[column];
2204 ind=UcolInd_+colBeg;
2205 indEnd=ind+UcolLengths_[column];
2206 uCol=Ucolumns_+colBeg;
2207 for ( ; ind!=indEnd; ++ind ){
2208 #if 0
2209 rhs1[ *ind ]-= (*uCol) * x1;
2210 rhs2[ *ind ]-= (*uCol) * x2;
2211 #else
2212 double value1=rhs1[ *ind ];
2213 rhs1[ *ind ] = value1 - (*uCol) * x1;
2214 double value2=rhs2[ *ind ];
2215 rhs2[ *ind ] = value2 - (*uCol) * x2;
2216 #endif
2217 ++uCol;
2218 }
2219 sol1[column]=x1;
2220 sol2[column]=x2;
2221 }
2222 }
2223 }
2224 for ( int k=numberSlacks_-1; k>=0; --k ){
2225 row=secRowOfU_[k];
2226 column=colOfU_[k];
2227 sol1[column]=-rhs1[row];
2228 sol2[column]=-rhs2[row];
2229 }
2230 }
2231
2232
xLeqb(double * b) const2233 void CoinSimpFactorization::xLeqb(double *b) const
2234 {
2235 double *rhs=b;
2236 int k, *ind, *indEnd, j;
2237 int colBeg;
2238 double x, *Lcol;
2239 // find last nonzero
2240 int last;
2241 for ( last=numberColumns_-1; last >= 0; --last ){
2242 if ( rhs[ rowOfU_[last] ] ) break;
2243 }
2244 // this seems to be faster
2245 if ( last >= 0 ){
2246 for ( j=last; j >=firstNumberSlacks_ ; --j ){
2247 k=rowOfU_[j];
2248 x=rhs[k];
2249 colBeg=LcolStarts_[k];
2250 ind=LcolInd_+colBeg;
2251 indEnd=ind+LcolLengths_[k];
2252 Lcol=Lcolumns_+colBeg;
2253 for ( ; ind!=indEnd; ++ind ){
2254 x -= (*Lcol) * rhs[ *ind ];
2255 ++Lcol;
2256 }
2257 rhs[k]=x;
2258 }
2259 } // if ( last >= 0 ){
2260 }
2261
2262
2263
xUeqb(double * b,double * sol) const2264 void CoinSimpFactorization::xUeqb(double *b, double *sol) const
2265 {
2266 double *rhs=b;
2267 int row, col, *ind, *indEnd, k;
2268 double xr;
2269 // now solve
2270 #if 1
2271 int rowBeg;
2272 double * uRow;
2273 for ( k=0; k<numberSlacks_; ++k ){
2274 row=secRowOfU_[k];
2275 col=colOfU_[k];
2276 xr=rhs[col];
2277 if ( xr!=0.0 ) {
2278 //if ( fabs(xr)> zeroTolerance_ ) {
2279 xr=-xr;
2280 rowBeg=UrowStarts_[row];
2281 ind=UrowInd_+rowBeg;
2282 indEnd=ind+UrowLengths_[row];
2283 uRow=Urows_+rowBeg;
2284 for ( ; ind!=indEnd; ++ind ){
2285 rhs[ *ind ]-= (*uRow) * xr;
2286 ++uRow;
2287 }
2288 sol[row]=xr;
2289 }
2290 else sol[row]=0.0;
2291 }
2292 for ( k=numberSlacks_; k<numberRows_; ++k ){
2293 row=secRowOfU_[k];
2294 col=colOfU_[k];
2295 xr=rhs[col];
2296 if ( xr!=0.0 ) {
2297 //if ( fabs(xr)> zeroTolerance_ ) {
2298 xr*=invOfPivots_[row];
2299 rowBeg=UrowStarts_[row];
2300 ind=UrowInd_+rowBeg;
2301 indEnd=ind+UrowLengths_[row];
2302 uRow=Urows_+rowBeg;
2303 for ( ; ind!=indEnd; ++ind ){
2304 rhs[ *ind ]-= (*uRow) * xr;
2305 ++uRow;
2306 }
2307 sol[row]=xr;
2308 }
2309 else sol[row]=0.0;
2310 }
2311 #else
2312 for ( k=0; k<numberSlacks_; ++k ){
2313 row=secRowOfU_[k];
2314 col=colOfU_[k];
2315 sol[row]=-rhs[col];
2316 }
2317 for ( k=numberSlacks_; k<numberRows_; ++k ){
2318 row=secRowOfU_[k];
2319 col=colOfU_[k];
2320 xr=rhs[col];
2321 int colBeg=UcolStarts_[col];
2322 ind=UcolInd_+colBeg;
2323 indEnd=ind+UcolLengths_[col];
2324 double * uCol=Ucolumns_+colBeg;
2325 for ( ; ind!=indEnd; ++ind,++uCol ){
2326 int iRow = *ind;
2327 double value = sol[iRow];
2328 double elementValue = *uCol;
2329 xr -= value*elementValue;
2330 }
2331 if ( xr!=0.0 ) {
2332 xr*=invOfPivots_[row];
2333 sol[row]=xr;
2334 }
2335 else sol[row]=0.0;
2336 }
2337 #endif
2338 }
2339
2340
2341
LUupdate(int newBasicCol)2342 int CoinSimpFactorization::LUupdate(int newBasicCol)
2343 {
2344 //checkU();
2345 // recover vector kept in ftran
2346 double *newColumn=vecKeep_;
2347 int *indNewColumn=indKeep_;
2348 int sizeNewColumn=keepSize_;
2349
2350 // remove elements of new column of U
2351 const int colBeg=UcolStarts_[newBasicCol];
2352 const int colEnd=colBeg+UcolLengths_[newBasicCol];
2353 for ( int i=colBeg; i<colEnd; ++i ){
2354 const int row=UcolInd_[i];
2355 const int colInRow=findInRow(row,newBasicCol);
2356 assert(colInRow >= 0);
2357 // remove from row
2358 const int rowEnd=UrowStarts_[row]+UrowLengths_[row];
2359 Urows_[colInRow]=Urows_[rowEnd-1];
2360 UrowInd_[colInRow]=UrowInd_[rowEnd-1];
2361 --UrowLengths_[row];
2362 }
2363 UcolLengths_[newBasicCol]=0;
2364 // now add new column to U
2365 int lastRowInU=-1;
2366 for ( int i=0; i < sizeNewColumn; ++i ){
2367 //if ( fabs(newColumn[i]) < zeroTolerance_ ) continue;
2368 const int row=indNewColumn[i];
2369 // add to row
2370 #ifdef COIN_SIMP_CAPACITY
2371 if ( UrowLengths_[row] + 1 > UrowCapacities_[row] )
2372 increaseRowSize(row, UrowLengths_[row] + 1);
2373 #endif
2374 const int rowEnd=UrowStarts_[row]+UrowLengths_[row];
2375 UrowInd_[rowEnd]=newBasicCol;
2376 Urows_[rowEnd]=newColumn[i];
2377 ++UrowLengths_[row];
2378 if ( lastRowInU < secRowPosition_[row] ) lastRowInU=secRowPosition_[row];
2379 }
2380 // add to Ucolumns
2381 #ifdef COIN_SIMP_CAPACITY
2382 if ( sizeNewColumn > UcolCapacities_[newBasicCol] )
2383 increaseColSize(newBasicCol, sizeNewColumn , true);
2384 #endif
2385 memcpy(&Ucolumns_[ UcolStarts_[newBasicCol] ], &newColumn[0],
2386 sizeNewColumn * sizeof(double) );
2387 memcpy(&UcolInd_[ UcolStarts_[newBasicCol] ], &indNewColumn[0],
2388 sizeNewColumn * sizeof(int) );
2389 UcolLengths_[newBasicCol]=sizeNewColumn;
2390
2391
2392 const int posNewCol=colPosition_[newBasicCol];
2393 if ( lastRowInU < posNewCol ){
2394 // matrix is singular
2395 return 1;
2396 }
2397 // permutations
2398 const int rowInU=secRowOfU_[posNewCol];
2399 const int colInU=colOfU_[posNewCol];
2400 for ( int i=posNewCol; i<lastRowInU; ++i ){
2401 int indx=secRowOfU_[i+1];
2402 secRowOfU_[i]=indx;
2403 secRowPosition_[indx]=i;
2404 int jndx=colOfU_[i+1];
2405 colOfU_[i]=jndx;
2406 colPosition_[jndx]=i;
2407 }
2408 secRowOfU_[lastRowInU]=rowInU;
2409 secRowPosition_[rowInU]=lastRowInU;
2410 colOfU_[lastRowInU]=colInU;
2411 colPosition_[colInU]=lastRowInU;
2412 if ( posNewCol < numberSlacks_ ){
2413 if ( lastRowInU >= numberSlacks_ )
2414 --numberSlacks_;
2415 else
2416 numberSlacks_= lastRowInU;
2417 }
2418 // rowInU will be transformed
2419 // denseVector_ is assumed to be initialized to zero
2420 const int rowBeg=UrowStarts_[rowInU];
2421 const int rowEnd=rowBeg+UrowLengths_[rowInU];
2422 for ( int i=rowBeg; i<rowEnd; ++i ){
2423 const int column=UrowInd_[i];
2424 denseVector_[column]=Urows_[i];
2425 // remove element
2426 const int indxRow=findInColumn(column,rowInU);
2427 assert( indxRow >= 0 );
2428 const int colEnd=UcolStarts_[column]+UcolLengths_[column];
2429 UcolInd_[indxRow]=UcolInd_[colEnd-1];
2430 Ucolumns_[indxRow]=Ucolumns_[colEnd-1];
2431 --UcolLengths_[column];
2432 }
2433 UrowLengths_[rowInU]=0;
2434 // rowInU is empty
2435 // increase Eta by (lastRowInU-posNewCol) elements
2436 newEta(rowInU, lastRowInU-posNewCol );
2437 assert(!EtaLengths_[lastEtaRow_]);
2438 int saveSize = EtaSize_;
2439 for ( int i=posNewCol; i<lastRowInU; ++i ){
2440 const int row=secRowOfU_[i];
2441 const int column=colOfU_[i];
2442 if ( denseVector_[column]==0.0 ) continue;
2443 const double multiplier=denseVector_[column]*invOfPivots_[row];
2444 denseVector_[column]=0.0;
2445 const int rowBeg=UrowStarts_[row];
2446 const int rowEnd=rowBeg+UrowLengths_[row];
2447 #if ARRAY
2448 for ( int j=rowBeg; j<rowEnd; ++j ){
2449 denseVector_[ UrowInd_[j] ]-= multiplier*Urows_[j];
2450 }
2451 #else
2452 int *ind=UrowInd_+rowBeg;
2453 int *indEnd=UrowInd_+rowEnd;
2454 double *uRow=Urows_+rowBeg;
2455 for ( ; ind!=indEnd; ++ind ){
2456 denseVector_[ *ind ]-= multiplier * (*uRow);
2457 ++uRow;
2458 }
2459 #endif
2460 // store multiplier
2461 Eta_[EtaSize_]=multiplier;
2462 EtaInd_[EtaSize_++]=row;
2463 }
2464 if (EtaSize_!=saveSize)
2465 EtaLengths_[lastEtaRow_]=EtaSize_ - saveSize;
2466 else
2467 --lastEtaRow_;
2468 // inverse of diagonal
2469 invOfPivots_[rowInU]=1.0/denseVector_[ colOfU_[lastRowInU] ];
2470 denseVector_[ colOfU_[lastRowInU] ]=0.0;
2471 // now store row
2472 int newEls=0;
2473 for ( int i=lastRowInU+1; i<numberColumns_; ++i ){
2474 const int column=colOfU_[i];
2475 const double coeff=denseVector_[column];
2476 denseVector_[column]=0.0;
2477 if ( fabs(coeff) < zeroTolerance_ ) continue;
2478 #ifdef COIN_SIMP_CAPACITY
2479 if ( UcolLengths_[column] + 1 > UcolCapacities_[column] ){
2480 increaseColSize(column, UcolLengths_[column] + 1, true);
2481 }
2482 #endif
2483 const int newInd=UcolStarts_[column]+UcolLengths_[column];
2484 UcolInd_[newInd]=rowInU;
2485 Ucolumns_[newInd]=coeff;
2486 ++UcolLengths_[column];
2487 workArea2_[newEls]=coeff;
2488 indVector_[newEls++]=column;
2489 }
2490 #ifdef COIN_SIMP_CAPACITY
2491 if ( UrowCapacities_[rowInU] < newEls )
2492 increaseRowSize(rowInU, newEls);
2493 #endif
2494 const int startRow=UrowStarts_[rowInU];
2495 memcpy(&Urows_[startRow],&workArea2_[0], newEls*sizeof(double) );
2496 memcpy(&UrowInd_[startRow],&indVector_[0], newEls*sizeof(int) );
2497 UrowLengths_[rowInU]=newEls;
2498 //
2499 if ( fabs( invOfPivots_[rowInU] ) > updateTol_ )
2500 return 2;
2501
2502 return 0;
2503 }
2504
2505
newEta(int row,int numNewElements)2506 void CoinSimpFactorization::newEta(int row, int numNewElements){
2507 if ( lastEtaRow_ == maxEtaRows_-1 ){
2508 int *iaux=new int[maxEtaRows_ + minIncrease_];
2509 memcpy(iaux, EtaPosition_, maxEtaRows_ * sizeof(int));
2510 delete [] EtaPosition_;
2511 EtaPosition_=iaux;
2512
2513 int *jaux=new int[maxEtaRows_ + minIncrease_];
2514 memcpy(jaux, EtaStarts_, maxEtaRows_ * sizeof(int));
2515 delete [] EtaStarts_;
2516 EtaStarts_=jaux;
2517
2518 int *kaux=new int[maxEtaRows_ + minIncrease_];
2519 memcpy(kaux, EtaLengths_, maxEtaRows_ * sizeof(int));
2520 delete [] EtaLengths_;
2521 EtaLengths_=kaux;
2522
2523 maxEtaRows_+=minIncrease_;
2524 }
2525 if ( EtaSize_ + numNewElements > EtaMaxCap_ ){
2526 int number= CoinMax(EtaSize_ + numNewElements - EtaMaxCap_, minIncrease_);
2527
2528 int *iaux=new int[EtaMaxCap_ + number];
2529 memcpy(iaux, EtaInd_, EtaSize_ * sizeof(int));
2530 delete [] EtaInd_;
2531 EtaInd_=iaux;
2532
2533 double *aux=new double[EtaMaxCap_ + number];
2534 memcpy(aux, Eta_, EtaSize_ * sizeof(double));
2535 delete [] Eta_;
2536 Eta_=aux;
2537
2538 EtaMaxCap_+=number;
2539 }
2540 EtaPosition_[++lastEtaRow_]=row;
2541 EtaStarts_[lastEtaRow_]=EtaSize_;
2542 EtaLengths_[lastEtaRow_]=0;
2543
2544 }
copyRowPermutations()2545 void CoinSimpFactorization::copyRowPermutations()
2546 {
2547 memcpy(&secRowOfU_[0], &rowOfU_[0],
2548 numberRows_ * sizeof(int) );
2549 memcpy(&secRowPosition_[0], &rowPosition_[0],
2550 numberRows_ * sizeof(int) );
2551 }
2552
Hxeqb(double * b) const2553 void CoinSimpFactorization::Hxeqb(double *b) const
2554 {
2555 double *rhs=b;
2556 int row, rowBeg, *ind, *indEnd;
2557 double xr, *eta;
2558 // now solve
2559 for ( int k=0; k <= lastEtaRow_; ++k ){
2560 row=EtaPosition_[k];
2561 rowBeg=EtaStarts_[k];
2562 xr=0.0;
2563 ind=EtaInd_+rowBeg;
2564 indEnd=ind+EtaLengths_[k];
2565 eta=Eta_+rowBeg;
2566 for ( ; ind!=indEnd; ++ind ){
2567 xr += rhs[ *ind ] * (*eta);
2568 ++eta;
2569 }
2570 rhs[row]-=xr;
2571 }
2572 }
2573
2574
2575
Hxeqb2(double * b1,double * b2) const2576 void CoinSimpFactorization::Hxeqb2(double *b1, double *b2) const
2577 {
2578 double *rhs1=b1;
2579 double *rhs2=b2;
2580 int row, rowBeg, *ind, *indEnd;
2581 double x1, x2, *eta;
2582 // now solve
2583 for ( int k=0; k <= lastEtaRow_; ++k ){
2584 row=EtaPosition_[k];
2585 rowBeg=EtaStarts_[k];
2586 x1=0.0;
2587 x2=0.0;
2588 ind=EtaInd_+rowBeg;
2589 indEnd=ind+EtaLengths_[k];
2590 eta=Eta_+rowBeg;
2591 for ( ; ind!=indEnd; ++ind ){
2592 x1 += rhs1[ *ind ] * (*eta);
2593 x2 += rhs2[ *ind ] * (*eta);
2594 ++eta;
2595 }
2596 rhs1[row]-=x1;
2597 rhs2[row]-=x2;
2598 }
2599 }
2600
2601
2602
2603
xHeqb(double * b) const2604 void CoinSimpFactorization::xHeqb(double *b) const
2605 {
2606 double *rhs=b;
2607 int row, rowBeg, *ind, *indEnd;
2608 double xr, *eta;
2609 // now solve
2610 for ( int k=lastEtaRow_; k >= 0; --k ){
2611 row=EtaPosition_[k];
2612 xr=rhs[row];
2613 if ( xr==0.0 ) continue;
2614 //if ( fabs(xr) <= zeroTolerance_ ) continue;
2615 rowBeg=EtaStarts_[k];
2616 ind=EtaInd_+rowBeg;
2617 indEnd=ind+EtaLengths_[k];
2618 eta=Eta_+rowBeg;
2619 for ( ; ind!=indEnd; ++ind ){
2620 rhs[ *ind ]-= xr * (*eta);
2621 ++eta;
2622 }
2623 }
2624 }
2625
2626
2627
ftran(double * b,double * sol,bool save) const2628 void CoinSimpFactorization::ftran(double *b, double *sol, bool save) const
2629 {
2630 Lxeqb(b);
2631 Hxeqb(b);
2632 if ( save ){
2633 // keep vector
2634 keepSize_=0;
2635 for ( int i=0; i<numberRows_; ++i ){
2636 if ( fabs(b[i]) < zeroTolerance_ ) continue;
2637 vecKeep_[keepSize_]=b[i];
2638 indKeep_[keepSize_++]=i;
2639 }
2640 }
2641 Uxeqb(b,sol);
2642 }
2643
ftran2(double * b1,double * sol1,double * b2,double * sol2) const2644 void CoinSimpFactorization::ftran2(double *b1, double *sol1, double *b2, double *sol2) const
2645 {
2646 Lxeqb2(b1,b2);
2647 Hxeqb2(b1,b2);
2648 // keep vector
2649 keepSize_=0;
2650 for ( int i=0; i<numberRows_; ++i ){
2651 if ( fabs(b1[i]) < zeroTolerance_ ) continue;
2652 vecKeep_[keepSize_]=b1[i];
2653 indKeep_[keepSize_++]=i;
2654 }
2655 Uxeqb2(b1,sol1,b2,sol2);
2656 }
2657
2658
2659
btran(double * b,double * sol) const2660 void CoinSimpFactorization::btran(double *b, double *sol) const
2661 {
2662 xUeqb(b, sol);
2663 xHeqb(sol);
2664 xLeqb(sol);
2665 }
2666