1 /* $Id: CoinAbcDenseFactorization.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2 // Copyright (C) 2008, International Business Machines
3 // Corporation and others, Copyright (C) 2012, FasterCoin. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include "CoinUtilsConfig.h"
7 #include "CoinPragma.hpp"
8
9 #include <cassert>
10 #include <cstdio>
11
12 #include "CoinAbcDenseFactorization.hpp"
13 #include "CoinAbcCommonFactorization.hpp"
14 #include "CoinIndexedVector.hpp"
15 #include "CoinHelperFunctions.hpp"
16 #include "CoinPackedMatrix.hpp"
17 #include "CoinFinite.hpp"
18 //:class CoinAbcDenseFactorization. Deals with Factorization and Updates
19 // CoinAbcDenseFactorization. Constructor
CoinAbcDenseFactorization()20 CoinAbcDenseFactorization::CoinAbcDenseFactorization()
21 : CoinAbcAnyFactorization()
22 {
23 gutsOfInitialize();
24 }
25
26 /// Copy constructor
CoinAbcDenseFactorization(const CoinAbcDenseFactorization & other)27 CoinAbcDenseFactorization::CoinAbcDenseFactorization(const CoinAbcDenseFactorization &other)
28 : CoinAbcAnyFactorization(other)
29 {
30 gutsOfInitialize();
31 gutsOfCopy(other);
32 }
33 // Clone
34 CoinAbcAnyFactorization *
clone() const35 CoinAbcDenseFactorization::clone() const
36 {
37 return new CoinAbcDenseFactorization(*this);
38 }
39 /// The real work of constructors etc
gutsOfDestructor()40 void CoinAbcDenseFactorization::gutsOfDestructor()
41 {
42 delete[] elements_;
43 delete[] pivotRow_;
44 delete[] workArea_;
45 elements_ = NULL;
46 pivotRow_ = NULL;
47 workArea_ = NULL;
48 numberRows_ = 0;
49 numberGoodU_ = 0;
50 status_ = -1;
51 maximumRows_ = 0;
52 #if ABC_PARALLEL == 2
53 parallelMode_ = 0;
54 #endif
55 maximumSpace_ = 0;
56 maximumRowsAdjusted_ = 0;
57 solveMode_ = 0;
58 }
gutsOfInitialize()59 void CoinAbcDenseFactorization::gutsOfInitialize()
60 {
61 pivotTolerance_ = 1.0e-1;
62 minimumPivotTolerance_ = 1.0e-1;
63 zeroTolerance_ = 1.0e-13;
64 areaFactor_ = 1.0;
65 numberDense_ = 0;
66 maximumPivots_ = 200;
67 relaxCheck_ = 1.0;
68 numberRows_ = 0;
69 numberGoodU_ = 0;
70 status_ = -1;
71 numberSlacks_ = 0;
72 numberPivots_ = 0;
73 maximumRows_ = 0;
74 #if ABC_PARALLEL == 2
75 parallelMode_ = 0;
76 #endif
77 maximumSpace_ = 0;
78 maximumRowsAdjusted_ = 0;
79 elements_ = NULL;
80 pivotRow_ = NULL;
81 workArea_ = NULL;
82 solveMode_ = 0;
83 }
84 // ~CoinAbcDenseFactorization. Destructor
~CoinAbcDenseFactorization()85 CoinAbcDenseFactorization::~CoinAbcDenseFactorization()
86 {
87 gutsOfDestructor();
88 }
89 // =
operator =(const CoinAbcDenseFactorization & other)90 CoinAbcDenseFactorization &CoinAbcDenseFactorization::operator=(const CoinAbcDenseFactorization &other)
91 {
92 if (this != &other) {
93 gutsOfDestructor();
94 gutsOfInitialize();
95 gutsOfCopy(other);
96 }
97 return *this;
98 }
99 #define WORK_MULT 2
gutsOfCopy(const CoinAbcDenseFactorization & other)100 void CoinAbcDenseFactorization::gutsOfCopy(const CoinAbcDenseFactorization &other)
101 {
102 pivotTolerance_ = other.pivotTolerance_;
103 minimumPivotTolerance_ = other.minimumPivotTolerance_;
104 zeroTolerance_ = other.zeroTolerance_;
105 areaFactor_ = other.areaFactor_;
106 numberDense_ = other.numberDense_;
107 relaxCheck_ = other.relaxCheck_;
108 numberRows_ = other.numberRows_;
109 maximumRows_ = other.maximumRows_;
110 #if ABC_PARALLEL == 2
111 parallelMode_ = other.parallelMode_;
112 #endif
113 maximumSpace_ = other.maximumSpace_;
114 maximumRowsAdjusted_ = other.maximumRowsAdjusted_;
115 solveMode_ = other.solveMode_;
116 numberGoodU_ = other.numberGoodU_;
117 maximumPivots_ = other.maximumPivots_;
118 numberPivots_ = other.numberPivots_;
119 factorElements_ = other.factorElements_;
120 status_ = other.status_;
121 numberSlacks_ = other.numberSlacks_;
122 if (other.pivotRow_) {
123 pivotRow_ = new int[2 * maximumRowsAdjusted_ + maximumPivots_];
124 CoinMemcpyN(other.pivotRow_, (2 * maximumRowsAdjusted_ + numberPivots_), pivotRow_);
125 elements_ = new CoinFactorizationDouble[maximumSpace_];
126 CoinMemcpyN(other.elements_, (maximumRowsAdjusted_ + numberPivots_) * maximumRowsAdjusted_, elements_);
127 workArea_ = new CoinFactorizationDouble[maximumRowsAdjusted_ * WORK_MULT];
128 CoinZeroN(workArea_, maximumRowsAdjusted_ * WORK_MULT);
129 } else {
130 elements_ = NULL;
131 pivotRow_ = NULL;
132 workArea_ = NULL;
133 }
134 }
135
136 // getAreas. Gets space for a factorization
137 //called by constructors
getAreas(int numberOfRows,int,CoinBigIndex,CoinBigIndex)138 void CoinAbcDenseFactorization::getAreas(int numberOfRows,
139 int /*numberOfColumns*/,
140 CoinBigIndex,
141 CoinBigIndex)
142 {
143
144 numberRows_ = numberOfRows;
145 numberDense_ = numberRows_;
146 if ((numberRows_ & (BLOCKING8 - 1)) != 0)
147 numberDense_ += (BLOCKING8 - (numberRows_ & (BLOCKING8 - 1)));
148 CoinBigIndex size = numberDense_ * (2 * numberDense_ + CoinMax(maximumPivots_, (numberDense_ + 1) >> 1));
149 if (size > maximumSpace_) {
150 delete[] elements_;
151 elements_ = new CoinFactorizationDouble[size];
152 maximumSpace_ = size;
153 }
154 if (numberRows_ > maximumRows_) {
155 maximumRows_ = numberRows_;
156 maximumRowsAdjusted_ = maximumRows_;
157 if ((maximumRows_ & (BLOCKING8 - 1)) != 0)
158 maximumRowsAdjusted_ += (BLOCKING8 - (maximumRows_ & (BLOCKING8 - 1)));
159 delete[] pivotRow_;
160 delete[] workArea_;
161 pivotRow_ = new int[2 * maximumRowsAdjusted_ + maximumPivots_];
162 workArea_ = new CoinFactorizationDouble[maximumRowsAdjusted_ * WORK_MULT];
163 }
164 }
165
166 // preProcess.
preProcess()167 void CoinAbcDenseFactorization::preProcess()
168 {
169 CoinBigIndex put = numberDense_ * numberDense_;
170 CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
171 CoinAbcMemset0(area, put);
172 int *indexRow = reinterpret_cast< int * >(elements_ + numberRows_ * numberRows_);
173 CoinBigIndex *starts = reinterpret_cast< CoinBigIndex * >(pivotRow_);
174 CoinFactorizationDouble *COIN_RESTRICT column = area;
175 //if (solveMode_==10)
176 //solveMode_=1;
177 if ((solveMode_ % 10) != 0) {
178 for (int i = 0; i < numberRows_; i++) {
179 for (CoinBigIndex j = starts[i]; j < starts[i + 1]; j++) {
180 int iRow = indexRow[j];
181 int iBlock = iRow >> 3;
182 iRow = iRow & (BLOCKING8 - 1);
183 column[iRow + BLOCKING8X8 * iBlock] = elements_[j];
184 }
185 if ((i & (BLOCKING8 - 1)) != (BLOCKING8 - 1))
186 column += BLOCKING8;
187 else
188 column += numberDense_ * BLOCKING8 - (BLOCKING8 - 1) * BLOCKING8;
189 }
190 for (int i = numberRows_; i < numberDense_; i++) {
191 int iRow = i;
192 int iBlock = iRow >> 3;
193 iRow = iRow & (BLOCKING8 - 1);
194 column[iRow + BLOCKING8X8 * iBlock] = 1.0;
195 if ((i & (BLOCKING8 - 1)) != (BLOCKING8 - 1))
196 column += BLOCKING8;
197 //else
198 //column += numberDense_*BLOCKING8-(BLOCKING8-1)*BLOCKING8;
199 }
200 } else {
201 // first or bad
202 for (int i = 0; i < numberRows_; i++) {
203 for (CoinBigIndex j = starts[i]; j < starts[i + 1]; j++) {
204 int iRow = indexRow[j];
205 column[iRow] = elements_[j];
206 }
207 column += numberDense_;
208 }
209 for (int i = numberRows_; i < numberDense_; i++) {
210 column[i] = 1.0;
211 column += numberDense_;
212 }
213 }
214 }
215
216 //Does factorization
factor(AbcSimplex *)217 int CoinAbcDenseFactorization::factor(AbcSimplex * /*model*/)
218 {
219 numberPivots_ = 0;
220 // ? take out
221 //printf("Debug non lapack dense factor\n");
222 //solveMode_&=~1;
223 CoinBigIndex put = numberDense_ * numberDense_;
224 CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
225 if ((solveMode_ % 10) != 0) {
226 // save last start
227 CoinBigIndex lastStart = pivotRow_[numberRows_];
228 status_ = CoinAbcDgetrf(numberDense_, numberDense_, area, numberDense_, pivotRow_ + numberDense_
229 #if ABC_PARALLEL == 2
230 ,
231 parallelMode_
232 #endif
233 );
234 // need to check size of pivots
235 if (!status_) {
236 // OK
237 solveMode_ = 1 + 10 * (solveMode_ / 10);
238 numberGoodU_ = numberRows_;
239 CoinZeroN(workArea_, 2 * numberRows_);
240 #if 0 //ndef NDEBUG
241 const CoinFactorizationDouble * column = area;
242 double smallest=COIN_DBL_MAX;
243 for (int i=0;i<numberRows_;i++) {
244 if (fabs(column[i])<smallest)
245 smallest = fabs(column[i]);
246 column += numberRows_;
247 }
248 if (smallest<1.0e-8)
249 printf("small el %g\n",smallest);
250 #endif
251 return 0;
252 } else {
253 solveMode_ = 10 * (solveMode_ / 10);
254 // redo matrix
255 // last start
256 pivotRow_[numberRows_] = lastStart;
257 preProcess();
258 }
259 }
260 status_ = 0;
261 for (int j = 0; j < numberRows_; j++) {
262 pivotRow_[j + numberDense_] = j;
263 }
264 CoinFactorizationDouble *elements = area;
265 numberGoodU_ = 0;
266 for (int i = 0; i < numberRows_; i++) {
267 int iRow = -1;
268 // Find largest
269 double largest = zeroTolerance_;
270 for (int j = i; j < numberRows_; j++) {
271 double value = fabs(elements[j]);
272 if (value > largest) {
273 largest = value;
274 iRow = j;
275 }
276 }
277 if (iRow >= 0) {
278 if (iRow != i) {
279 // swap
280 assert(iRow > i);
281 CoinFactorizationDouble *elementsA = area;
282 for (int k = 0; k <= i; k++) {
283 // swap
284 CoinFactorizationDouble value = elementsA[i];
285 elementsA[i] = elementsA[iRow];
286 elementsA[iRow] = value;
287 elementsA += numberDense_;
288 }
289 int iPivot = pivotRow_[i + numberDense_];
290 pivotRow_[i + numberDense_] = pivotRow_[iRow + numberDense_];
291 pivotRow_[iRow + numberDense_] = iPivot;
292 }
293 CoinFactorizationDouble pivotValue = 1.0 / elements[i];
294 elements[i] = pivotValue;
295 for (int j = i + 1; j < numberRows_; j++) {
296 elements[j] *= pivotValue;
297 }
298 // Update rest
299 CoinFactorizationDouble *elementsA = elements;
300 for (int k = i + 1; k < numberRows_; k++) {
301 elementsA += numberDense_;
302 // swap
303 if (iRow != i) {
304 CoinFactorizationDouble value = elementsA[i];
305 elementsA[i] = elementsA[iRow];
306 elementsA[iRow] = value;
307 }
308 CoinFactorizationDouble value = elementsA[i];
309 for (int j = i + 1; j < numberRows_; j++) {
310 elementsA[j] -= value * elements[j];
311 }
312 }
313 } else {
314 status_ = -1;
315 break;
316 }
317 numberGoodU_++;
318 elements += numberDense_;
319 }
320 for (int j = 0; j < numberRows_; j++) {
321 int k = pivotRow_[j + numberDense_];
322 pivotRow_[k] = j;
323 }
324 //assert (status_);
325 //if (!status_)
326 //solveMode_=10; // for next time
327 //for (int j=0;j<numberRows_;j++) {
328 //int k = pivotRow_[j+numberDense_];
329 //pivotRow_[k]=j;
330 //}
331 return status_;
332 }
333 // Makes a non-singular basis by replacing variables
makeNonSingular(int * sequence)334 void CoinAbcDenseFactorization::makeNonSingular(int *sequence)
335 {
336 // Replace bad ones by correct slack
337 int *workArea = reinterpret_cast< int * >(workArea_);
338 int i;
339 for (i = 0; i < numberRows_; i++)
340 workArea[i] = -1;
341 for (i = 0; i < numberGoodU_; i++) {
342 int iOriginal = pivotRow_[i + numberDense_];
343 workArea[iOriginal] = i;
344 //workArea[i]=iOriginal;
345 }
346 int lastRow = -1;
347 for (i = 0; i < numberRows_; i++) {
348 if (workArea[i] == -1) {
349 lastRow = i;
350 break;
351 }
352 }
353 assert(lastRow >= 0);
354 for (i = numberGoodU_; i < numberRows_; i++) {
355 assert(lastRow < numberRows_);
356 // Put slack in basis
357 sequence[i] = lastRow;
358 lastRow++;
359 for (; lastRow < numberRows_; lastRow++) {
360 if (workArea[lastRow] == -1)
361 break;
362 }
363 }
364 }
365 //#define DENSE_PERMUTE
366 // Does post processing on valid factorization - putting variables on correct rows
postProcess(const int * sequence,int * pivotVariable)367 void CoinAbcDenseFactorization::postProcess(const int *sequence, int *pivotVariable)
368 {
369 if ((solveMode_ % 10) == 0) {
370 for (int i = 0; i < numberRows_; i++) {
371 int k = sequence[i];
372 #ifdef DENSE_PERMUTE
373 pivotVariable[pivotRow_[i + numberDense_]] = k;
374 #else
375 //pivotVariable[pivotRow_[i]]=k;
376 //pivotVariable[pivotRow_[i]]=k;
377 pivotVariable[i] = k;
378 k = pivotRow_[i];
379 pivotRow_[i] = pivotRow_[i + numberDense_];
380 pivotRow_[i + numberDense_] = k;
381 #endif
382 }
383 } else {
384 // lapack
385 for (int i = 0; i < numberRows_; i++) {
386 int k = sequence[i];
387 pivotVariable[i] = k;
388 }
389 }
390 }
391 /* Replaces one Column to basis,
392 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
393 partial update already in U */
replaceColumn(CoinIndexedVector * regionSparse,int pivotRow,double pivotCheck,bool,double)394 int CoinAbcDenseFactorization::replaceColumn(CoinIndexedVector *regionSparse,
395 int pivotRow,
396 double pivotCheck,
397 bool /*skipBtranU*/,
398 double /*acceptablePivot*/)
399 {
400 if (numberPivots_ == maximumPivots_)
401 return 3;
402 CoinFactorizationDouble *elements = elements_ + numberDense_ * numberPivots_;
403 double *region = regionSparse->denseVector();
404 int *regionIndex = regionSparse->getIndices();
405 int numberNonZero = regionSparse->getNumElements();
406 int i;
407 memset(elements, 0, numberRows_ * sizeof(CoinFactorizationDouble));
408 CoinFactorizationDouble pivotValue = pivotCheck;
409 if (fabs(pivotValue) < zeroTolerance_)
410 return 2;
411 pivotValue = 1.0 / pivotValue;
412 if ((solveMode_ % 10) == 0) {
413 if (regionSparse->packedMode()) {
414 for (i = 0; i < numberNonZero; i++) {
415 int iRow = regionIndex[i];
416 double value = region[i];
417 #ifdef DENSE_PERMUTE
418 iRow = pivotRow_[iRow]; // permute
419 #endif
420 elements[iRow] = value;
421 ;
422 }
423 } else {
424 // not packed! - from user pivot?
425 for (i = 0; i < numberNonZero; i++) {
426 int iRow = regionIndex[i];
427 double value = region[iRow];
428 #ifdef DENSE_PERMUTE
429 iRow = pivotRow_[iRow]; // permute
430 #endif
431 elements[iRow] = value;
432 ;
433 }
434 }
435 int realPivotRow = pivotRow_[pivotRow];
436 elements[realPivotRow] = pivotValue;
437 pivotRow_[2 * numberDense_ + numberPivots_] = realPivotRow;
438 } else {
439 // lapack
440 if (regionSparse->packedMode()) {
441 for (i = 0; i < numberNonZero; i++) {
442 int iRow = regionIndex[i];
443 double value = region[i];
444 elements[iRow] = value;
445 ;
446 }
447 } else {
448 // not packed! - from user pivot?
449 for (i = 0; i < numberNonZero; i++) {
450 int iRow = regionIndex[i];
451 double value = region[iRow];
452 elements[iRow] = value;
453 ;
454 }
455 }
456 elements[pivotRow] = pivotValue;
457 pivotRow_[2 * numberDense_ + numberPivots_] = pivotRow;
458 }
459 numberPivots_++;
460 return 0;
461 }
462 /* Checks if can replace one Column to basis,
463 returns 0=OK, 1=Probably OK, 2=singular, 3=no room, 5 max pivots */
checkReplacePart2(int pivotRow,double btranAlpha,double ftranAlpha,long double ftAlpha,double acceptablePivot)464 int CoinAbcDenseFactorization::checkReplacePart2(int pivotRow,
465 double btranAlpha,
466 double ftranAlpha,
467 #ifdef ABC_LONG_FACTORIZATION
468 long
469 #endif
470 double ftAlpha,
471 double acceptablePivot)
472 {
473 if (numberPivots_ == maximumPivots_)
474 return 3;
475 if (fabs(ftranAlpha) < zeroTolerance_)
476 return 2;
477 return 0;
478 }
479 /* Replaces one Column to basis,
480 partial update already in U */
replaceColumnPart3(const AbcSimplex * model,CoinIndexedVector * regionSparse,CoinIndexedVector * tableauColumn,int pivotRow,long double alpha)481 void CoinAbcDenseFactorization::replaceColumnPart3(const AbcSimplex *model,
482 CoinIndexedVector *regionSparse,
483 CoinIndexedVector *tableauColumn,
484 int pivotRow,
485 #ifdef ABC_LONG_FACTORIZATION
486 long
487 #endif
488 double alpha)
489 {
490 CoinFactorizationDouble *elements = elements_ + numberDense_ * numberPivots_;
491 double *region = tableauColumn->denseVector();
492 int *regionIndex = tableauColumn->getIndices();
493 int numberNonZero = tableauColumn->getNumElements();
494 int i;
495 memset(elements, 0, numberRows_ * sizeof(CoinFactorizationDouble));
496 double pivotValue = 1.0 / alpha;
497 if ((solveMode_ % 10) == 0) {
498 if (tableauColumn->packedMode()) {
499 for (i = 0; i < numberNonZero; i++) {
500 int iRow = regionIndex[i];
501 double value = region[i];
502 #ifdef DENSE_PERMUTE
503 iRow = pivotRow_[iRow]; // permute
504 #endif
505 elements[iRow] = value;
506 ;
507 }
508 } else {
509 // not packed! - from user pivot?
510 for (i = 0; i < numberNonZero; i++) {
511 int iRow = regionIndex[i];
512 double value = region[iRow];
513 #ifdef DENSE_PERMUTE
514 iRow = pivotRow_[iRow]; // permute
515 #endif
516 elements[iRow] = value;
517 ;
518 }
519 }
520 int realPivotRow = pivotRow_[pivotRow];
521 elements[realPivotRow] = pivotValue;
522 pivotRow_[2 * numberDense_ + numberPivots_] = realPivotRow;
523 } else {
524 // lapack
525 if (tableauColumn->packedMode()) {
526 for (i = 0; i < numberNonZero; i++) {
527 int iRow = regionIndex[i];
528 double value = region[i];
529 elements[iRow] = value;
530 ;
531 }
532 } else {
533 // not packed! - from user pivot?
534 for (i = 0; i < numberNonZero; i++) {
535 int iRow = regionIndex[i];
536 double value = region[iRow];
537 elements[iRow] = value;
538 ;
539 }
540 }
541 elements[pivotRow] = pivotValue;
542 pivotRow_[2 * numberDense_ + numberPivots_] = pivotRow;
543 }
544 numberPivots_++;
545 }
546 static void
CoinAbcDlaswp1(double * COIN_RESTRICT a,int m,int * ipiv)547 CoinAbcDlaswp1(double *COIN_RESTRICT a, int m, int *ipiv)
548 {
549 for (int i = 0; i < m; i++) {
550 int ip = ipiv[i];
551 if (ip != i) {
552 double temp = a[i];
553 a[i] = a[ip];
554 a[ip] = temp;
555 }
556 }
557 }
558 static void
CoinAbcDlaswp1Backwards(double * COIN_RESTRICT a,int m,int * ipiv)559 CoinAbcDlaswp1Backwards(double *COIN_RESTRICT a, int m, int *ipiv)
560 {
561 for (int i = m - 1; i >= 0; i--) {
562 int ip = ipiv[i];
563 if (ip != i) {
564 double temp = a[i];
565 a[i] = a[ip];
566 a[ip] = temp;
567 }
568 }
569 }
570 /* This version has same effect as above with FTUpdate==false
571 so number returned is always >=0 */
updateColumn(CoinIndexedVector & regionSparse) const572 int CoinAbcDenseFactorization::updateColumn(CoinIndexedVector ®ionSparse) const
573 {
574 double *region = regionSparse.denseVector();
575 CoinBigIndex put = numberDense_ * numberDense_;
576 CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
577 CoinFactorizationDouble *COIN_RESTRICT elements = area;
578 if ((solveMode_ % 10) == 0) {
579 // base factorization L
580 for (int i = 0; i < numberRows_; i++) {
581 double value = region[i];
582 for (int j = i + 1; j < numberRows_; j++) {
583 region[j] -= value * elements[j];
584 }
585 elements += numberDense_;
586 }
587 elements = area + numberRows_ * numberDense_;
588 // base factorization U
589 for (int i = numberRows_ - 1; i >= 0; i--) {
590 elements -= numberDense_;
591 CoinFactorizationDouble value = region[i] * elements[i];
592 region[i] = value;
593 for (int j = 0; j < i; j++) {
594 region[j] -= value * elements[j];
595 }
596 }
597 } else {
598 CoinAbcDlaswp1(region, numberRows_, pivotRow_ + numberDense_);
599 CoinAbcDgetrs('N', numberDense_, area, region);
600 }
601 // now updates
602 elements = elements_;
603 for (int i = 0; i < numberPivots_; i++) {
604 int iPivot = pivotRow_[i + 2 * numberDense_];
605 CoinFactorizationDouble value = region[iPivot] * elements[iPivot];
606 for (int j = 0; j < numberRows_; j++) {
607 region[j] -= value * elements[j];
608 }
609 region[iPivot] = value;
610 elements += numberDense_;
611 }
612 regionSparse.setNumElements(0);
613 regionSparse.scan(0, numberRows_);
614 return 0;
615 }
616
617 /* Updates one column for dual steepest edge weights (FTRAN) */
updateWeights(CoinIndexedVector & regionSparse) const618 void CoinAbcDenseFactorization::updateWeights(CoinIndexedVector ®ionSparse) const
619 {
620 updateColumn(regionSparse);
621 }
622
updateTwoColumnsFT(CoinIndexedVector & regionSparse2,CoinIndexedVector & regionSparse3)623 int CoinAbcDenseFactorization::updateTwoColumnsFT(CoinIndexedVector ®ionSparse2,
624 CoinIndexedVector ®ionSparse3)
625 {
626 updateColumn(regionSparse2);
627 updateColumn(regionSparse3);
628 return 0;
629 }
630
631 /* Updates one column (BTRAN) from regionSparse2
632 regionSparse starts as zero and is zero at end
633 Note - if regionSparse2 packed on input - will be packed on output
634 */
updateColumnTranspose(CoinIndexedVector & regionSparse2) const635 int CoinAbcDenseFactorization::updateColumnTranspose(CoinIndexedVector ®ionSparse2) const
636 {
637 double *region = regionSparse2.denseVector();
638 CoinFactorizationDouble *elements = elements_ + numberDense_ * numberPivots_;
639 // updates
640 for (int i = numberPivots_ - 1; i >= 0; i--) {
641 elements -= numberDense_;
642 int iPivot = pivotRow_[i + 2 * numberDense_];
643 CoinFactorizationDouble value = region[iPivot]; //*elements[iPivot];
644 for (int j = 0; j < iPivot; j++) {
645 value -= region[j] * elements[j];
646 }
647 for (int j = iPivot + 1; j < numberRows_; j++) {
648 value -= region[j] * elements[j];
649 }
650 region[iPivot] = value * elements[iPivot];
651 }
652 CoinBigIndex put = numberDense_ * numberDense_;
653 CoinFactorizationDouble *COIN_RESTRICT area = elements_ + maximumSpace_ - put;
654 if ((solveMode_ % 10) == 0) {
655 // base factorization U
656 elements = area;
657 for (int i = 0; i < numberRows_; i++) {
658 //CoinFactorizationDouble value = region[i]*elements[i];
659 CoinFactorizationDouble value = region[i];
660 for (int j = 0; j < i; j++) {
661 value -= region[j] * elements[j];
662 }
663 //region[i] = value;
664 region[i] = value * elements[i];
665 elements += numberDense_;
666 }
667 // base factorization L
668 elements = area + numberDense_ * numberRows_;
669 for (int i = numberRows_ - 1; i >= 0; i--) {
670 elements -= numberDense_;
671 CoinFactorizationDouble value = region[i];
672 for (int j = i + 1; j < numberRows_; j++) {
673 value -= region[j] * elements[j];
674 }
675 region[i] = value;
676 }
677 } else {
678 CoinAbcDgetrs('T', numberDense_, area, region);
679 CoinAbcDlaswp1Backwards(region, numberRows_, pivotRow_ + numberDense_);
680 }
681 regionSparse2.setNumElements(0);
682 regionSparse2.scan(0, numberRows_);
683 return 0;
684 }
685 /* Updates one column (FTRAN) */
updateColumnCpu(CoinIndexedVector & regionSparse,int whichCpu) const686 void CoinAbcAnyFactorization::updateColumnCpu(CoinIndexedVector ®ionSparse, int whichCpu) const
687 {
688 updateColumn(regionSparse);
689 }
690 /* Updates one column (BTRAN) */
updateColumnTransposeCpu(CoinIndexedVector & regionSparse,int whichCpu) const691 void CoinAbcAnyFactorization::updateColumnTransposeCpu(CoinIndexedVector ®ionSparse, int whichCpu) const
692 {
693 updateColumnTranspose(regionSparse);
694 }
695 // Default constructor
CoinAbcAnyFactorization()696 CoinAbcAnyFactorization::CoinAbcAnyFactorization()
697 : pivotTolerance_(1.0e-1)
698 , minimumPivotTolerance_(1.0e-1)
699 , zeroTolerance_(1.0e-13)
700 , relaxCheck_(1.0)
701 , factorElements_(0)
702 , numberRows_(0)
703 , numberGoodU_(0)
704 , maximumPivots_(200)
705 , numberPivots_(0)
706 , status_(-1)
707 , solveMode_(0)
708 {
709 }
710 // Copy constructor
CoinAbcAnyFactorization(const CoinAbcAnyFactorization & other)711 CoinAbcAnyFactorization::CoinAbcAnyFactorization(const CoinAbcAnyFactorization &other)
712 : pivotTolerance_(other.pivotTolerance_)
713 , minimumPivotTolerance_(other.minimumPivotTolerance_)
714 , zeroTolerance_(other.zeroTolerance_)
715 , relaxCheck_(other.relaxCheck_)
716 , factorElements_(other.factorElements_)
717 , numberRows_(other.numberRows_)
718 , numberGoodU_(other.numberGoodU_)
719 , maximumPivots_(other.maximumPivots_)
720 , numberPivots_(other.numberPivots_)
721 , status_(other.status_)
722 , solveMode_(other.solveMode_)
723 {
724 }
725 // Destructor
~CoinAbcAnyFactorization()726 CoinAbcAnyFactorization::~CoinAbcAnyFactorization()
727 {
728 }
729 // = copy
operator =(const CoinAbcAnyFactorization & other)730 CoinAbcAnyFactorization &CoinAbcAnyFactorization::operator=(const CoinAbcAnyFactorization &other)
731 {
732 if (this != &other) {
733 pivotTolerance_ = other.pivotTolerance_;
734 minimumPivotTolerance_ = other.minimumPivotTolerance_;
735 zeroTolerance_ = other.zeroTolerance_;
736 relaxCheck_ = other.relaxCheck_;
737 factorElements_ = other.factorElements_;
738 numberRows_ = other.numberRows_;
739 numberGoodU_ = other.numberGoodU_;
740 maximumPivots_ = other.maximumPivots_;
741 numberPivots_ = other.numberPivots_;
742 status_ = other.status_;
743 solveMode_ = other.solveMode_;
744 }
745 return *this;
746 }
pivotTolerance(double value)747 void CoinAbcAnyFactorization::pivotTolerance(double value)
748 {
749 if (value > 0.0 && value <= 1.0) {
750 pivotTolerance_ = CoinMax(value, minimumPivotTolerance_);
751 }
752 }
minimumPivotTolerance(double value)753 void CoinAbcAnyFactorization::minimumPivotTolerance(double value)
754 {
755 if (value > 0.0 && value <= 1.0) {
756 minimumPivotTolerance_ = value;
757 }
758 }
zeroTolerance(double value)759 void CoinAbcAnyFactorization::zeroTolerance(double value)
760 {
761 if (value > 0.0 && value < 1.0) {
762 zeroTolerance_ = value;
763 }
764 }
maximumPivots(int value)765 void CoinAbcAnyFactorization::maximumPivots(int value)
766 {
767 if (value > maximumPivots_) {
768 delete[] pivotRow_;
769 int n = maximumRows_;
770 if ((n & (BLOCKING8 - 1)) != 0)
771 n += (BLOCKING8 - (n & (BLOCKING8 - 1)));
772 pivotRow_ = new int[2 * n + value];
773 }
774 maximumPivots_ = value;
775 }
776 // Number of entries in each row
numberInRow() const777 int *CoinAbcAnyFactorization::numberInRow() const
778 {
779 return reinterpret_cast< int * >(workArea_);
780 }
781 // Number of entries in each column
numberInColumn() const782 int *CoinAbcAnyFactorization::numberInColumn() const
783 {
784 return (reinterpret_cast< int * >(workArea_)) + numberRows_;
785 }
786 // Returns array to put basis starts in
787 CoinBigIndex *
starts() const788 CoinAbcAnyFactorization::starts() const
789 {
790 return reinterpret_cast< CoinBigIndex * >(pivotRow_);
791 }
792 // Returns array to put basis elements in
793 CoinFactorizationDouble *
elements() const794 CoinAbcAnyFactorization::elements() const
795 {
796 return elements_;
797 }
798 // Returns pivot row
pivotRow() const799 int *CoinAbcAnyFactorization::pivotRow() const
800 {
801 return pivotRow_;
802 }
803 // Returns work area
804 CoinFactorizationDouble *
workArea() const805 CoinAbcAnyFactorization::workArea() const
806 {
807 return workArea_;
808 }
809 // Returns int work area
intWorkArea() const810 int *CoinAbcAnyFactorization::intWorkArea() const
811 {
812 return reinterpret_cast< int * >(workArea_);
813 }
814 // Returns permute back
permuteBack() const815 int *CoinAbcAnyFactorization::permuteBack() const
816 {
817 return pivotRow_ + numberRows_;
818 }
819 // Returns pivot column
pivotColumn() const820 int *CoinAbcAnyFactorization::pivotColumn() const
821 {
822 return permute();
823 }
824 // Returns true if wants tableauColumn in replaceColumn
wantsTableauColumn() const825 bool CoinAbcAnyFactorization::wantsTableauColumn() const
826 {
827 return true;
828 }
829 /* Useful information for factorization
830 0 - iteration number
831 whereFrom is 0 for factorize and 1 for replaceColumn
832 */
setUsefulInformation(const int *,int)833 void CoinAbcAnyFactorization::setUsefulInformation(const int *, int)
834 {
835 }
836
837 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
838 */
839