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