1 /* $Id: CoinDenseFactorization.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
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 #include "CoinPragma.hpp"
8
9 #include <cassert>
10 #include <cstdio>
11
12 #include "CoinDenseFactorization.hpp"
13 #include "CoinIndexedVector.hpp"
14 #include "CoinHelperFunctions.hpp"
15 #include "CoinPackedMatrix.hpp"
16 #include "CoinFinite.hpp"
17 #if COIN_BIG_DOUBLE == 1
18 #undef COIN_FACTORIZATION_DENSE_CODE
19 #endif
20 #ifdef COIN_FACTORIZATION_DENSE_CODE
21 // using simple lapack interface
22 extern "C" {
23 /** LAPACK Fortran subroutine DGETRF. */
24 void F77_FUNC(dgetrf, DGETRF)(ipfint *m, ipfint *n,
25 double *A, ipfint *ldA,
26 ipfint *ipiv, ipfint *info);
27 /** LAPACK Fortran subroutine DGETRS. */
28 void F77_FUNC(dgetrs, DGETRS)(char *trans, cipfint *n,
29 cipfint *nrhs, const double *A, cipfint *ldA,
30 cipfint *ipiv, double *B, cipfint *ldB, ipfint *info,
31 int trans_len);
32 }
33 #endif
34 //:class CoinDenseFactorization. Deals with Factorization and Updates
35 // CoinDenseFactorization. Constructor
CoinDenseFactorization()36 CoinDenseFactorization::CoinDenseFactorization()
37 : CoinOtherFactorization()
38 {
39 gutsOfInitialize();
40 }
41
42 /// Copy constructor
CoinDenseFactorization(const CoinDenseFactorization & other)43 CoinDenseFactorization::CoinDenseFactorization(const CoinDenseFactorization &other)
44 : CoinOtherFactorization(other)
45 {
46 gutsOfInitialize();
47 gutsOfCopy(other);
48 }
49 // Clone
50 CoinOtherFactorization *
clone() const51 CoinDenseFactorization::clone() const
52 {
53 return new CoinDenseFactorization(*this);
54 }
55 /// The real work of constructors etc
gutsOfDestructor()56 void CoinDenseFactorization::gutsOfDestructor()
57 {
58 delete[] elements_;
59 delete[] pivotRow_;
60 delete[] workArea_;
61 elements_ = NULL;
62 pivotRow_ = NULL;
63 workArea_ = NULL;
64 numberRows_ = 0;
65 numberColumns_ = 0;
66 numberGoodU_ = 0;
67 status_ = -1;
68 maximumRows_ = 0;
69 maximumSpace_ = 0;
70 solveMode_ = 0;
71 }
gutsOfInitialize()72 void CoinDenseFactorization::gutsOfInitialize()
73 {
74 pivotTolerance_ = 1.0e-1;
75 zeroTolerance_ = 1.0e-13;
76 #ifndef COIN_FAST_CODE
77 slackValue_ = -1.0;
78 #endif
79 maximumPivots_ = 200;
80 relaxCheck_ = 1.0;
81 numberRows_ = 0;
82 numberColumns_ = 0;
83 numberGoodU_ = 0;
84 status_ = -1;
85 numberPivots_ = 0;
86 maximumRows_ = 0;
87 maximumSpace_ = 0;
88 elements_ = NULL;
89 pivotRow_ = NULL;
90 workArea_ = NULL;
91 solveMode_ = 0;
92 }
93 // ~CoinDenseFactorization. Destructor
~CoinDenseFactorization()94 CoinDenseFactorization::~CoinDenseFactorization()
95 {
96 gutsOfDestructor();
97 }
98 // =
operator =(const CoinDenseFactorization & other)99 CoinDenseFactorization &CoinDenseFactorization::operator=(const CoinDenseFactorization &other)
100 {
101 if (this != &other) {
102 gutsOfDestructor();
103 gutsOfInitialize();
104 gutsOfCopy(other);
105 }
106 return *this;
107 }
108 #ifdef COIN_FACTORIZATION_DENSE_CODE
109 #define WORK_MULT 2
110 #else
111 #define WORK_MULT 2
112 #endif
gutsOfCopy(const CoinDenseFactorization & other)113 void CoinDenseFactorization::gutsOfCopy(const CoinDenseFactorization &other)
114 {
115 pivotTolerance_ = other.pivotTolerance_;
116 zeroTolerance_ = other.zeroTolerance_;
117 #ifndef COIN_FAST_CODE
118 slackValue_ = other.slackValue_;
119 #endif
120 relaxCheck_ = other.relaxCheck_;
121 numberRows_ = other.numberRows_;
122 numberColumns_ = other.numberColumns_;
123 maximumRows_ = other.maximumRows_;
124 maximumSpace_ = other.maximumSpace_;
125 solveMode_ = other.solveMode_;
126 numberGoodU_ = other.numberGoodU_;
127 maximumPivots_ = other.maximumPivots_;
128 numberPivots_ = other.numberPivots_;
129 factorElements_ = other.factorElements_;
130 status_ = other.status_;
131 if (other.pivotRow_) {
132 pivotRow_ = new int[2 * maximumRows_ + maximumPivots_];
133 CoinMemcpyN(other.pivotRow_, (2 * maximumRows_ + numberPivots_), pivotRow_);
134 elements_ = new CoinFactorizationDouble[maximumSpace_];
135 CoinMemcpyN(other.elements_, (maximumRows_ + numberPivots_) * maximumRows_, elements_);
136 workArea_ = new CoinFactorizationDouble[maximumRows_ * WORK_MULT];
137 CoinZeroN(workArea_, maximumRows_ * WORK_MULT);
138 } else {
139 elements_ = NULL;
140 pivotRow_ = NULL;
141 workArea_ = NULL;
142 }
143 }
144
145 // getAreas. Gets space for a factorization
146 //called by constructors
getAreas(int numberOfRows,int numberOfColumns,int,int)147 void CoinDenseFactorization::getAreas(int numberOfRows,
148 int numberOfColumns,
149 int,
150 int)
151 {
152
153 numberRows_ = numberOfRows;
154 numberColumns_ = numberOfColumns;
155 int size = numberRows_ * (numberRows_ + CoinMax(maximumPivots_, (numberRows_ + 1) >> 1));
156 if (size > maximumSpace_) {
157 delete[] elements_;
158 elements_ = new CoinFactorizationDouble[size];
159 maximumSpace_ = size;
160 }
161 if (numberRows_ > maximumRows_) {
162 maximumRows_ = numberRows_;
163 delete[] pivotRow_;
164 delete[] workArea_;
165 pivotRow_ = new int[2 * maximumRows_ + maximumPivots_];
166 workArea_ = new CoinFactorizationDouble[maximumRows_ * WORK_MULT];
167 }
168 }
169
170 // preProcess.
preProcess()171 void CoinDenseFactorization::preProcess()
172 {
173 // could do better than this but this only a demo
174 int put = numberRows_ * numberRows_;
175 int *indexRow = reinterpret_cast< int * >(elements_ + put);
176 int *starts = reinterpret_cast< int * >(pivotRow_);
177 put = numberRows_ * numberColumns_;
178 for (int i = numberColumns_ - 1; i >= 0; i--) {
179 put -= numberRows_;
180 memset(workArea_, 0, numberRows_ * sizeof(CoinFactorizationDouble));
181 assert(starts[i] <= put);
182 for (int j = starts[i]; j < starts[i + 1]; j++) {
183 int iRow = indexRow[j];
184 workArea_[iRow] = elements_[j];
185 }
186 // move to correct position
187 CoinMemcpyN(workArea_, numberRows_, elements_ + put);
188 }
189 }
190
191 //Does factorization
factor()192 int CoinDenseFactorization::factor()
193 {
194 numberPivots_ = 0;
195 status_ = 0;
196 #ifdef COIN_FACTORIZATION_DENSE_CODE
197 if (numberRows_ == numberColumns_ && (solveMode_ % 10) != 0) {
198 int info;
199 F77_FUNC(dgetrf, DGETRF)
200 (&numberRows_, &numberRows_,
201 elements_, &numberRows_, pivotRow_,
202 &info);
203 // need to check size of pivots
204 if (!info) {
205 // OK
206 solveMode_ = 1 + 10 * (solveMode_ / 10);
207 numberGoodU_ = numberRows_;
208 CoinZeroN(workArea_, 2 * numberRows_);
209 #if 0 //ndef NDEBUG
210 const CoinFactorizationDouble * column = elements_;
211 double smallest=COIN_DBL_MAX;
212 for (int i=0;i<numberRows_;i++) {
213 if (fabs(column[i])<smallest)
214 smallest = fabs(column[i]);
215 column += numberRows_;
216 }
217 if (smallest<1.0e-8)
218 printf("small el %g\n",smallest);
219 #endif
220 return 0;
221 } else {
222 solveMode_ = 10 * (solveMode_ / 10);
223 }
224 }
225 #endif
226 for (int j = 0; j < numberRows_; j++) {
227 pivotRow_[j + numberRows_] = j;
228 }
229 CoinFactorizationDouble *elements = elements_;
230 numberGoodU_ = 0;
231 for (int i = 0; i < numberColumns_; i++) {
232 int iRow = -1;
233 // Find largest
234 double largest = zeroTolerance_;
235 for (int j = i; j < numberRows_; j++) {
236 double value = fabs(elements[j]);
237 if (value > largest) {
238 largest = value;
239 iRow = j;
240 }
241 }
242 if (iRow >= 0) {
243 if (iRow != i) {
244 // swap
245 assert(iRow > i);
246 CoinFactorizationDouble *elementsA = elements_;
247 for (int k = 0; k <= i; k++) {
248 // swap
249 CoinFactorizationDouble value = elementsA[i];
250 elementsA[i] = elementsA[iRow];
251 elementsA[iRow] = value;
252 elementsA += numberRows_;
253 }
254 int iPivot = pivotRow_[i + numberRows_];
255 pivotRow_[i + numberRows_] = pivotRow_[iRow + numberRows_];
256 pivotRow_[iRow + numberRows_] = iPivot;
257 }
258 CoinFactorizationDouble pivotValue = 1.0 / elements[i];
259 elements[i] = pivotValue;
260 for (int j = i + 1; j < numberRows_; j++) {
261 elements[j] *= pivotValue;
262 }
263 // Update rest
264 CoinFactorizationDouble *elementsA = elements;
265 for (int k = i + 1; k < numberColumns_; k++) {
266 elementsA += numberRows_;
267 // swap
268 if (iRow != i) {
269 CoinFactorizationDouble value = elementsA[i];
270 elementsA[i] = elementsA[iRow];
271 elementsA[iRow] = value;
272 }
273 CoinFactorizationDouble value = elementsA[i];
274 for (int j = i + 1; j < numberRows_; j++) {
275 elementsA[j] -= value * elements[j];
276 }
277 }
278 } else {
279 status_ = -1;
280 break;
281 }
282 numberGoodU_++;
283 elements += numberRows_;
284 }
285 for (int j = 0; j < numberRows_; j++) {
286 int k = pivotRow_[j + numberRows_];
287 pivotRow_[k] = j;
288 }
289 return status_;
290 }
291 // Makes a non-singular basis by replacing variables
makeNonSingular(int * sequence,int numberColumns)292 void CoinDenseFactorization::makeNonSingular(int *sequence, int numberColumns)
293 {
294 // Replace bad ones by correct slack
295 int *workArea = reinterpret_cast< int * >(workArea_);
296 int i;
297 for (i = 0; i < numberRows_; i++)
298 workArea[i] = -1;
299 for (i = 0; i < numberGoodU_; i++) {
300 int iOriginal = pivotRow_[i + numberRows_];
301 workArea[iOriginal] = i;
302 //workArea[i]=iOriginal;
303 }
304 int lastRow = -1;
305 for (i = 0; i < numberRows_; i++) {
306 if (workArea[i] == -1) {
307 lastRow = i;
308 break;
309 }
310 }
311 assert(lastRow >= 0);
312 for (i = numberGoodU_; i < numberRows_; i++) {
313 assert(lastRow < numberRows_);
314 // Put slack in basis
315 sequence[i] = lastRow + numberColumns;
316 lastRow++;
317 for (; lastRow < numberRows_; lastRow++) {
318 if (workArea[lastRow] == -1)
319 break;
320 }
321 }
322 }
323 #define DENSE_PERMUTE
324 // Does post processing on valid factorization - putting variables on correct rows
postProcess(const int * sequence,int * pivotVariable)325 void CoinDenseFactorization::postProcess(const int *sequence, int *pivotVariable)
326 {
327 #ifdef COIN_FACTORIZATION_DENSE_CODE
328 if ((solveMode_ % 10) == 0) {
329 #endif
330 for (int i = 0; i < numberRows_; i++) {
331 int k = sequence[i];
332 #ifdef DENSE_PERMUTE
333 pivotVariable[pivotRow_[i + numberRows_]] = k;
334 #else
335 //pivotVariable[pivotRow_[i]]=k;
336 //pivotVariable[pivotRow_[i]]=k;
337 pivotVariable[i] = k;
338 k = pivotRow_[i];
339 pivotRow_[i] = pivotRow_[i + numberRows_];
340 pivotRow_[i + numberRows_] = k;
341 #endif
342 }
343 #ifdef COIN_FACTORIZATION_DENSE_CODE
344 } else {
345 // lapack
346 for (int i = 0; i < numberRows_; i++) {
347 int k = sequence[i];
348 pivotVariable[i] = k;
349 }
350 }
351 #endif
352 }
353 /* Replaces one Column to basis,
354 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
355 If checkBeforeModifying is true will do all accuracy checks
356 before modifying factorization. Whether to set this depends on
357 speed considerations. You could just do this on first iteration
358 after factorization and thereafter re-factorize
359 partial update already in U */
replaceColumn(CoinIndexedVector * regionSparse,int pivotRow,double pivotCheck,bool,double)360 int CoinDenseFactorization::replaceColumn(CoinIndexedVector *regionSparse,
361 int pivotRow,
362 double pivotCheck,
363 bool /*checkBeforeModifying*/,
364 double /*acceptablePivot*/)
365 {
366 if (numberPivots_ == maximumPivots_)
367 return 3;
368 CoinFactorizationDouble *elements = elements_ + numberRows_ * (numberColumns_ + numberPivots_);
369 double *region = regionSparse->denseVector();
370 int *regionIndex = regionSparse->getIndices();
371 int numberNonZero = regionSparse->getNumElements();
372 int i;
373 memset(elements, 0, numberRows_ * sizeof(CoinFactorizationDouble));
374 CoinFactorizationDouble pivotValue = pivotCheck;
375 if (fabs(pivotValue) < zeroTolerance_)
376 return 2;
377 pivotValue = 1.0 / pivotValue;
378 #ifdef COIN_FACTORIZATION_DENSE_CODE
379 if ((solveMode_ % 10) == 0) {
380 #endif
381 if (regionSparse->packedMode()) {
382 for (i = 0; i < numberNonZero; i++) {
383 int iRow = regionIndex[i];
384 double value = region[i];
385 #ifdef DENSE_PERMUTE
386 iRow = pivotRow_[iRow]; // permute
387 #endif
388 elements[iRow] = value;
389 ;
390 }
391 } else {
392 // not packed! - from user pivot?
393 for (i = 0; i < numberNonZero; i++) {
394 int iRow = regionIndex[i];
395 double value = region[iRow];
396 #ifdef DENSE_PERMUTE
397 iRow = pivotRow_[iRow]; // permute
398 #endif
399 elements[iRow] = value;
400 ;
401 }
402 }
403 int realPivotRow = pivotRow_[pivotRow];
404 elements[realPivotRow] = pivotValue;
405 pivotRow_[2 * numberRows_ + numberPivots_] = realPivotRow;
406 #ifdef COIN_FACTORIZATION_DENSE_CODE
407 } else {
408 // lapack
409 if (regionSparse->packedMode()) {
410 for (i = 0; i < numberNonZero; i++) {
411 int iRow = regionIndex[i];
412 double value = region[i];
413 elements[iRow] = value;
414 ;
415 }
416 } else {
417 // not packed! - from user pivot?
418 for (i = 0; i < numberNonZero; i++) {
419 int iRow = regionIndex[i];
420 double value = region[iRow];
421 elements[iRow] = value;
422 ;
423 }
424 }
425 elements[pivotRow] = pivotValue;
426 pivotRow_[2 * numberRows_ + numberPivots_] = pivotRow;
427 }
428 #endif
429 numberPivots_++;
430 return 0;
431 }
432 /* This version has same effect as above with FTUpdate==false
433 so number returned is always >=0 */
updateColumn(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2,bool noPermute) const434 int CoinDenseFactorization::updateColumn(CoinIndexedVector *regionSparse,
435 CoinIndexedVector *regionSparse2,
436 bool noPermute) const
437 {
438 assert(numberRows_ == numberColumns_);
439 double *region2 = regionSparse2->denseVector();
440 int *regionIndex = regionSparse2->getIndices();
441 int numberNonZero = regionSparse2->getNumElements();
442 double *region = regionSparse->denseVector();
443 #ifdef COIN_FACTORIZATION_DENSE_CODE
444 if ((solveMode_ % 10) == 0) {
445 #endif
446 if (!regionSparse2->packedMode()) {
447 if (!noPermute) {
448 for (int j = 0; j < numberRows_; j++) {
449 int iRow = pivotRow_[j + numberRows_];
450 region[j] = region2[iRow];
451 region2[iRow] = 0.0;
452 }
453 } else {
454 // can't due to check mode assert (regionSparse==regionSparse2);
455 region = regionSparse2->denseVector();
456 }
457 } else {
458 // packed mode
459 assert(!noPermute);
460 for (int j = 0; j < numberNonZero; j++) {
461 int jRow = regionIndex[j];
462 int iRow = pivotRow_[jRow];
463 region[iRow] = region2[j];
464 region2[j] = 0.0;
465 }
466 }
467 #ifdef COIN_FACTORIZATION_DENSE_CODE
468 } else {
469 // lapack
470 if (!regionSparse2->packedMode()) {
471 if (!noPermute) {
472 for (int j = 0; j < numberRows_; j++) {
473 region[j] = region2[j];
474 region2[j] = 0.0;
475 }
476 } else {
477 // can't due to check mode assert (regionSparse==regionSparse2);
478 region = regionSparse2->denseVector();
479 }
480 } else {
481 // packed mode
482 assert(!noPermute);
483 for (int j = 0; j < numberNonZero; j++) {
484 int jRow = regionIndex[j];
485 region[jRow] = region2[j];
486 region2[j] = 0.0;
487 }
488 }
489 }
490 #endif
491 int i;
492 CoinFactorizationDouble *elements = elements_;
493 #ifdef COIN_FACTORIZATION_DENSE_CODE
494 if ((solveMode_ % 10) == 0) {
495 #endif
496 // base factorization L
497 for (i = 0; i < numberColumns_; i++) {
498 double value = region[i];
499 for (int j = i + 1; j < numberRows_; j++) {
500 region[j] -= value * elements[j];
501 }
502 elements += numberRows_;
503 }
504 elements = elements_ + numberRows_ * numberRows_;
505 // base factorization U
506 for (i = numberColumns_ - 1; i >= 0; i--) {
507 elements -= numberRows_;
508 CoinFactorizationDouble value = region[i] * elements[i];
509 region[i] = value;
510 for (int j = 0; j < i; j++) {
511 region[j] -= value * elements[j];
512 }
513 }
514 #ifdef COIN_FACTORIZATION_DENSE_CODE
515 } else {
516 char trans = 'N';
517 int ione = 1;
518 int info;
519 F77_FUNC(dgetrs, DGETRS)
520 (&trans, &numberRows_, &ione, elements_, &numberRows_,
521 pivotRow_, region, &numberRows_, &info, 1);
522 }
523 #endif
524 // now updates
525 elements = elements_ + numberRows_ * numberRows_;
526 for (i = 0; i < numberPivots_; i++) {
527 int iPivot = pivotRow_[i + 2 * numberRows_];
528 CoinFactorizationDouble value = region[iPivot] * elements[iPivot];
529 for (int j = 0; j < numberRows_; j++) {
530 region[j] -= value * elements[j];
531 }
532 region[iPivot] = value;
533 elements += numberRows_;
534 }
535 // permute back and get nonzeros
536 numberNonZero = 0;
537 #ifdef COIN_FACTORIZATION_DENSE_CODE
538 if ((solveMode_ % 10) == 0) {
539 #endif
540 if (!noPermute) {
541 if (!regionSparse2->packedMode()) {
542 for (int j = 0; j < numberRows_; j++) {
543 #ifdef DENSE_PERMUTE
544 int iRow = pivotRow_[j];
545 #else
546 int iRow = j;
547 #endif
548 double value = region[iRow];
549 region[iRow] = 0.0;
550 if (fabs(value) > zeroTolerance_) {
551 region2[j] = value;
552 regionIndex[numberNonZero++] = j;
553 }
554 }
555 } else {
556 // packed mode
557 for (int j = 0; j < numberRows_; j++) {
558 #ifdef DENSE_PERMUTE
559 int iRow = pivotRow_[j];
560 #else
561 int iRow = j;
562 #endif
563 double value = region[iRow];
564 region[iRow] = 0.0;
565 if (fabs(value) > zeroTolerance_) {
566 region2[numberNonZero] = value;
567 regionIndex[numberNonZero++] = j;
568 }
569 }
570 }
571 } else {
572 for (int j = 0; j < numberRows_; j++) {
573 double value = region[j];
574 if (fabs(value) > zeroTolerance_) {
575 regionIndex[numberNonZero++] = j;
576 } else {
577 region[j] = 0.0;
578 }
579 }
580 }
581 #ifdef COIN_FACTORIZATION_DENSE_CODE
582 } else {
583 // lapack
584 if (!noPermute) {
585 if (!regionSparse2->packedMode()) {
586 for (int j = 0; j < numberRows_; j++) {
587 double value = region[j];
588 region[j] = 0.0;
589 if (fabs(value) > zeroTolerance_) {
590 region2[j] = value;
591 regionIndex[numberNonZero++] = j;
592 }
593 }
594 } else {
595 // packed mode
596 for (int j = 0; j < numberRows_; j++) {
597 double value = region[j];
598 region[j] = 0.0;
599 if (fabs(value) > zeroTolerance_) {
600 region2[numberNonZero] = value;
601 regionIndex[numberNonZero++] = j;
602 }
603 }
604 }
605 } else {
606 for (int j = 0; j < numberRows_; j++) {
607 double value = region[j];
608 if (fabs(value) > zeroTolerance_) {
609 regionIndex[numberNonZero++] = j;
610 } else {
611 region[j] = 0.0;
612 }
613 }
614 }
615 }
616 #endif
617 regionSparse2->setNumElements(numberNonZero);
618 return 0;
619 }
620
updateTwoColumnsFT(CoinIndexedVector * regionSparse1,CoinIndexedVector * regionSparse2,CoinIndexedVector * regionSparse3,bool)621 int CoinDenseFactorization::updateTwoColumnsFT(CoinIndexedVector *regionSparse1,
622 CoinIndexedVector *regionSparse2,
623 CoinIndexedVector *regionSparse3,
624 bool /*noPermute*/)
625 {
626 #ifdef COIN_FACTORIZATION_DENSE_CODE
627 #if 0
628 CoinIndexedVector s2(*regionSparse2);
629 CoinIndexedVector s3(*regionSparse3);
630 updateColumn(regionSparse1,&s2);
631 updateColumn(regionSparse1,&s3);
632 #endif
633 if ((solveMode_ % 10) == 0) {
634 #endif
635 updateColumn(regionSparse1, regionSparse2);
636 updateColumn(regionSparse1, regionSparse3);
637 #ifdef COIN_FACTORIZATION_DENSE_CODE
638 } else {
639 // lapack
640 assert(numberRows_ == numberColumns_);
641 double *region2 = regionSparse2->denseVector();
642 int *regionIndex2 = regionSparse2->getIndices();
643 int numberNonZero2 = regionSparse2->getNumElements();
644 CoinFactorizationDouble *regionW2 = workArea_;
645 if (!regionSparse2->packedMode()) {
646 for (int j = 0; j < numberRows_; j++) {
647 regionW2[j] = region2[j];
648 region2[j] = 0.0;
649 }
650 } else {
651 // packed mode
652 for (int j = 0; j < numberNonZero2; j++) {
653 int jRow = regionIndex2[j];
654 regionW2[jRow] = region2[j];
655 region2[j] = 0.0;
656 }
657 }
658 double *region3 = regionSparse3->denseVector();
659 int *regionIndex3 = regionSparse3->getIndices();
660 int numberNonZero3 = regionSparse3->getNumElements();
661 CoinFactorizationDouble *regionW3 = workArea_ + numberRows_;
662 if (!regionSparse3->packedMode()) {
663 for (int j = 0; j < numberRows_; j++) {
664 regionW3[j] = region3[j];
665 region3[j] = 0.0;
666 }
667 } else {
668 // packed mode
669 for (int j = 0; j < numberNonZero3; j++) {
670 int jRow = regionIndex3[j];
671 regionW3[jRow] = region3[j];
672 region3[j] = 0.0;
673 }
674 }
675 int i;
676 CoinFactorizationDouble *elements = elements_;
677 char trans = 'N';
678 int itwo = 2;
679 int info;
680 F77_FUNC(dgetrs, DGETRS)
681 (&trans, &numberRows_, &itwo, elements_, &numberRows_,
682 pivotRow_, workArea_, &numberRows_, &info, 1);
683 // now updates
684 elements = elements_ + numberRows_ * numberRows_;
685 for (i = 0; i < numberPivots_; i++) {
686 int iPivot = pivotRow_[i + 2 * numberRows_];
687 CoinFactorizationDouble value2 = regionW2[iPivot] * elements[iPivot];
688 CoinFactorizationDouble value3 = regionW3[iPivot] * elements[iPivot];
689 for (int j = 0; j < numberRows_; j++) {
690 regionW2[j] -= value2 * elements[j];
691 regionW3[j] -= value3 * elements[j];
692 }
693 regionW2[iPivot] = value2;
694 regionW3[iPivot] = value3;
695 elements += numberRows_;
696 }
697 // permute back and get nonzeros
698 numberNonZero2 = 0;
699 if (!regionSparse2->packedMode()) {
700 for (int j = 0; j < numberRows_; j++) {
701 double value = regionW2[j];
702 regionW2[j] = 0.0;
703 if (fabs(value) > zeroTolerance_) {
704 region2[j] = value;
705 regionIndex2[numberNonZero2++] = j;
706 }
707 }
708 } else {
709 // packed mode
710 for (int j = 0; j < numberRows_; j++) {
711 double value = regionW2[j];
712 regionW2[j] = 0.0;
713 if (fabs(value) > zeroTolerance_) {
714 region2[numberNonZero2] = value;
715 regionIndex2[numberNonZero2++] = j;
716 }
717 }
718 }
719 regionSparse2->setNumElements(numberNonZero2);
720 numberNonZero3 = 0;
721 if (!regionSparse3->packedMode()) {
722 for (int j = 0; j < numberRows_; j++) {
723 double value = regionW3[j];
724 regionW3[j] = 0.0;
725 if (fabs(value) > zeroTolerance_) {
726 region3[j] = value;
727 regionIndex3[numberNonZero3++] = j;
728 }
729 }
730 } else {
731 // packed mode
732 for (int j = 0; j < numberRows_; j++) {
733 double value = regionW3[j];
734 regionW3[j] = 0.0;
735 if (fabs(value) > zeroTolerance_) {
736 region3[numberNonZero3] = value;
737 regionIndex3[numberNonZero3++] = j;
738 }
739 }
740 }
741 regionSparse3->setNumElements(numberNonZero3);
742 #if 0
743 printf("Good2==\n");
744 s2.print();
745 printf("Bad2==\n");
746 regionSparse2->print();
747 printf("======\n");
748 printf("Good3==\n");
749 s3.print();
750 printf("Bad3==\n");
751 regionSparse3->print();
752 printf("======\n");
753 #endif
754 }
755 #endif
756 return 0;
757 }
758
759 /* Updates one column (BTRAN) from regionSparse2
760 regionSparse starts as zero and is zero at end
761 Note - if regionSparse2 packed on input - will be packed on output
762 */
updateColumnTranspose(CoinIndexedVector * regionSparse,CoinIndexedVector * regionSparse2) const763 int CoinDenseFactorization::updateColumnTranspose(CoinIndexedVector *regionSparse,
764 CoinIndexedVector *regionSparse2) const
765 {
766 assert(numberRows_ == numberColumns_);
767 double *region2 = regionSparse2->denseVector();
768 int *regionIndex = regionSparse2->getIndices();
769 int numberNonZero = regionSparse2->getNumElements();
770 double *region = regionSparse->denseVector();
771 #ifdef COIN_FACTORIZATION_DENSE_CODE
772 if ((solveMode_ % 10) == 0) {
773 #endif
774 if (!regionSparse2->packedMode()) {
775 for (int j = 0; j < numberRows_; j++) {
776 #ifdef DENSE_PERMUTE
777 int iRow = pivotRow_[j];
778 #else
779 int iRow = j;
780 #endif
781 region[iRow] = region2[j];
782 region2[j] = 0.0;
783 }
784 } else {
785 for (int j = 0; j < numberNonZero; j++) {
786 int jRow = regionIndex[j];
787 #ifdef DENSE_PERMUTE
788 int iRow = pivotRow_[jRow];
789 #else
790 int iRow = jRow;
791 #endif
792 region[iRow] = region2[j];
793 region2[j] = 0.0;
794 }
795 }
796 #ifdef COIN_FACTORIZATION_DENSE_CODE
797 } else {
798 // lapack
799 if (!regionSparse2->packedMode()) {
800 for (int j = 0; j < numberRows_; j++) {
801 region[j] = region2[j];
802 region2[j] = 0.0;
803 }
804 } else {
805 for (int j = 0; j < numberNonZero; j++) {
806 int jRow = regionIndex[j];
807 region[jRow] = region2[j];
808 region2[j] = 0.0;
809 }
810 }
811 }
812 #endif
813 int i;
814 CoinFactorizationDouble *elements = elements_ + numberRows_ * (numberRows_ + numberPivots_);
815 // updates
816 for (i = numberPivots_ - 1; i >= 0; i--) {
817 elements -= numberRows_;
818 int iPivot = pivotRow_[i + 2 * numberRows_];
819 CoinFactorizationDouble value = region[iPivot]; //*elements[iPivot];
820 for (int j = 0; j < iPivot; j++) {
821 value -= region[j] * elements[j];
822 }
823 for (int j = iPivot + 1; j < numberRows_; j++) {
824 value -= region[j] * elements[j];
825 }
826 region[iPivot] = value * elements[iPivot];
827 }
828 #ifdef COIN_FACTORIZATION_DENSE_CODE
829 if ((solveMode_ % 10) == 0) {
830 #endif
831 // base factorization U
832 elements = elements_;
833 for (i = 0; i < numberColumns_; i++) {
834 //CoinFactorizationDouble value = region[i]*elements[i];
835 CoinFactorizationDouble value = region[i];
836 for (int j = 0; j < i; j++) {
837 value -= region[j] * elements[j];
838 }
839 //region[i] = value;
840 region[i] = value * elements[i];
841 elements += numberRows_;
842 }
843 // base factorization L
844 elements = elements_ + numberRows_ * numberRows_;
845 for (i = numberColumns_ - 1; i >= 0; i--) {
846 elements -= numberRows_;
847 CoinFactorizationDouble value = region[i];
848 for (int j = i + 1; j < numberRows_; j++) {
849 value -= region[j] * elements[j];
850 }
851 region[i] = value;
852 }
853 #ifdef COIN_FACTORIZATION_DENSE_CODE
854 } else {
855 char trans = 'T';
856 int ione = 1;
857 int info;
858 F77_FUNC(dgetrs, DGETRS)
859 (&trans, &numberRows_, &ione, elements_, &numberRows_,
860 pivotRow_, region, &numberRows_, &info, 1);
861 }
862 #endif
863 // permute back and get nonzeros
864 numberNonZero = 0;
865 #ifdef COIN_FACTORIZATION_DENSE_CODE
866 if ((solveMode_ % 10) == 0) {
867 #endif
868 if (!regionSparse2->packedMode()) {
869 for (int j = 0; j < numberRows_; j++) {
870 int iRow = pivotRow_[j + numberRows_];
871 double value = region[j];
872 region[j] = 0.0;
873 if (fabs(value) > zeroTolerance_) {
874 region2[iRow] = value;
875 regionIndex[numberNonZero++] = iRow;
876 }
877 }
878 } else {
879 for (int j = 0; j < numberRows_; j++) {
880 int iRow = pivotRow_[j + numberRows_];
881 double value = region[j];
882 region[j] = 0.0;
883 if (fabs(value) > zeroTolerance_) {
884 region2[numberNonZero] = value;
885 regionIndex[numberNonZero++] = iRow;
886 }
887 }
888 }
889 #ifdef COIN_FACTORIZATION_DENSE_CODE
890 } else {
891 // lapack
892 if (!regionSparse2->packedMode()) {
893 for (int j = 0; j < numberRows_; j++) {
894 double value = region[j];
895 region[j] = 0.0;
896 if (fabs(value) > zeroTolerance_) {
897 region2[j] = value;
898 regionIndex[numberNonZero++] = j;
899 }
900 }
901 } else {
902 for (int j = 0; j < numberRows_; j++) {
903 double value = region[j];
904 region[j] = 0.0;
905 if (fabs(value) > zeroTolerance_) {
906 region2[numberNonZero] = value;
907 regionIndex[numberNonZero++] = j;
908 }
909 }
910 }
911 }
912 #endif
913 regionSparse2->setNumElements(numberNonZero);
914 return 0;
915 }
916 // Default constructor
CoinOtherFactorization()917 CoinOtherFactorization::CoinOtherFactorization()
918 : pivotTolerance_(1.0e-1)
919 , zeroTolerance_(1.0e-13)
920 ,
921 #ifndef COIN_FAST_CODE
922 slackValue_(-1.0)
923 ,
924 #endif
925 relaxCheck_(1.0)
926 , factorElements_(0)
927 , numberRows_(0)
928 , numberColumns_(0)
929 , numberGoodU_(0)
930 , maximumPivots_(200)
931 , numberPivots_(0)
932 , status_(-1)
933 , solveMode_(0)
934 {
935 }
936 // Copy constructor
CoinOtherFactorization(const CoinOtherFactorization & other)937 CoinOtherFactorization::CoinOtherFactorization(const CoinOtherFactorization &other)
938 : pivotTolerance_(other.pivotTolerance_)
939 , zeroTolerance_(other.zeroTolerance_)
940 ,
941 #ifndef COIN_FAST_CODE
942 slackValue_(other.slackValue_)
943 ,
944 #endif
945 relaxCheck_(other.relaxCheck_)
946 , factorElements_(other.factorElements_)
947 , numberRows_(other.numberRows_)
948 , numberColumns_(other.numberColumns_)
949 , numberGoodU_(other.numberGoodU_)
950 , maximumPivots_(other.maximumPivots_)
951 , numberPivots_(other.numberPivots_)
952 , status_(other.status_)
953 , solveMode_(other.solveMode_)
954 {
955 }
956 // Destructor
~CoinOtherFactorization()957 CoinOtherFactorization::~CoinOtherFactorization()
958 {
959 }
960 // = copy
operator =(const CoinOtherFactorization & other)961 CoinOtherFactorization &CoinOtherFactorization::operator=(const CoinOtherFactorization &other)
962 {
963 if (this != &other) {
964 pivotTolerance_ = other.pivotTolerance_;
965 zeroTolerance_ = other.zeroTolerance_;
966 #ifndef COIN_FAST_CODE
967 slackValue_ = other.slackValue_;
968 #endif
969 relaxCheck_ = other.relaxCheck_;
970 factorElements_ = other.factorElements_;
971 numberRows_ = other.numberRows_;
972 numberColumns_ = other.numberColumns_;
973 numberGoodU_ = other.numberGoodU_;
974 maximumPivots_ = other.maximumPivots_;
975 numberPivots_ = other.numberPivots_;
976 status_ = other.status_;
977 solveMode_ = other.solveMode_;
978 }
979 return *this;
980 }
pivotTolerance(double value)981 void CoinOtherFactorization::pivotTolerance(double value)
982 {
983 if (value > 0.0 && value <= 1.0) {
984 pivotTolerance_ = value;
985 }
986 }
zeroTolerance(double value)987 void CoinOtherFactorization::zeroTolerance(double value)
988 {
989 if (value > 0.0 && value < 1.0) {
990 zeroTolerance_ = value;
991 }
992 }
993 #ifndef COIN_FAST_CODE
slackValue(double value)994 void CoinOtherFactorization::slackValue(double value)
995 {
996 if (value >= 0.0) {
997 slackValue_ = 1.0;
998 } else {
999 slackValue_ = -1.0;
1000 }
1001 }
1002 #endif
maximumPivots(int value)1003 void CoinOtherFactorization::maximumPivots(int value)
1004 {
1005 if (value > maximumPivots_) {
1006 delete[] pivotRow_;
1007 pivotRow_ = new int[2 * maximumRows_ + value];
1008 }
1009 maximumPivots_ = value;
1010 }
1011 // Number of entries in each row
numberInRow() const1012 int *CoinOtherFactorization::numberInRow() const
1013 {
1014 return reinterpret_cast< int * >(workArea_);
1015 }
1016 // Number of entries in each column
numberInColumn() const1017 int *CoinOtherFactorization::numberInColumn() const
1018 {
1019 return (reinterpret_cast< int * >(workArea_)) + numberRows_;
1020 }
1021 // Returns array to put basis starts in
starts() const1022 int *CoinOtherFactorization::starts() const
1023 {
1024 return reinterpret_cast< int * >(pivotRow_);
1025 }
1026 // Returns array to put basis elements in
1027 CoinFactorizationDouble *
elements() const1028 CoinOtherFactorization::elements() const
1029 {
1030 return elements_;
1031 }
1032 // Returns pivot row
pivotRow() const1033 int *CoinOtherFactorization::pivotRow() const
1034 {
1035 return pivotRow_;
1036 }
1037 // Returns work area
1038 CoinFactorizationDouble *
workArea() const1039 CoinOtherFactorization::workArea() const
1040 {
1041 return workArea_;
1042 }
1043 // Returns int work area
intWorkArea() const1044 int *CoinOtherFactorization::intWorkArea() const
1045 {
1046 return reinterpret_cast< int * >(workArea_);
1047 }
1048 // Returns permute back
permuteBack() const1049 int *CoinOtherFactorization::permuteBack() const
1050 {
1051 return pivotRow_ + numberRows_;
1052 }
1053 // Returns true if wants tableauColumn in replaceColumn
wantsTableauColumn() const1054 bool CoinOtherFactorization::wantsTableauColumn() const
1055 {
1056 return true;
1057 }
1058 /* Useful information for factorization
1059 0 - iteration number
1060 whereFrom is 0 for factorize and 1 for replaceColumn
1061 */
setUsefulInformation(const int *,int)1062 void CoinOtherFactorization::setUsefulInformation(const int *, int)
1063 {
1064 }
1065
1066 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1067 */
1068