CoinUtils 2.11.4
Loading...
Searching...
No Matches
CoinFactorization.hpp
Go to the documentation of this file.
1/* $Id$ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6/*
7 Authors
8
9 John Forrest
10
11 */
12#ifndef CoinFactorization_H
13#define CoinFactorization_H
14#define EXTRA_U_SPACE 4
15//#define COIN_ONE_ETA_COPY 100
16
17#include <iostream>
18#include <string>
19#include <cassert>
20#include <cstdio>
21#include <cmath>
22#include "CoinTypes.hpp"
23#include "CoinIndexedVector.hpp"
24
51 friend void CoinFactorizationUnitTest(const std::string &mpsDir);
52
53public:
60
66 void show_self() const;
68 int saveFactorization(const char *file) const;
72 int restoreFactorization(const char *file, bool factor = false);
74 void sort() const;
78
88 int factorize(const CoinPackedMatrix &matrix,
89 int rowIsBasic[], int columnIsBasic[],
90 double areaFactor = 0.0);
102 int numberColumns,
103 int numberElements,
104 int maximumL,
105 int maximumU,
106 const int indicesRow[],
107 const int indicesColumn[], const double elements[],
108 int permutation[],
109 double areaFactor = 0.0);
115 int numberColumns,
116 int estimateNumberElements,
117 int *indicesRow[],
118 int *indicesColumn[],
119 CoinFactorizationDouble *elements[],
120 double areaFactor = 0.0);
127 int factorizePart2(int permutation[], int exactNumberElements);
129 double conditionNumber() const;
130
132
136 inline int status() const
137 {
138 return status_;
139 }
141 inline void setStatus(int value)
142 {
143 status_ = value;
144 }
146 inline int pivots() const
147 {
148 return numberPivots_;
149 }
151 inline void setPivots(int value)
152 {
153 numberPivots_ = value;
154 }
156 inline int *permute() const
157 {
158 return permute_.array();
159 }
161 inline int *pivotColumn() const
162 {
163 return pivotColumn_.array();
164 }
167 {
168 return pivotRegion_.array();
169 }
171 inline int *permuteBack() const
172 {
173 return permuteBack_.array();
174 }
176 inline int *lastRow() const
177 {
178 return lastRow_.array();
179 }
182 inline int *pivotColumnBack() const
183 {
184 //return firstCount_.array();
185 return pivotColumnBack_.array();
186 }
188 inline int *startRowL() const
189 {
190 return startRowL_.array();
191 }
192
194 inline int *startColumnL() const
195 {
196 return startColumnL_.array();
197 }
198
200 inline int *indexColumnL() const
201 {
202 return indexColumnL_.array();
203 }
204
206 inline int *indexRowL() const
207 {
208 return indexRowL_.array();
209 }
210
213 {
214 return elementByRowL_.array();
215 }
216
218 inline int numberRowsExtra() const
219 {
220 return numberRowsExtra_;
221 }
223 inline void setNumberRows(int value)
224 {
225 numberRows_ = value;
226 }
228 inline int numberRows() const
229 {
230 return numberRows_;
231 }
233 inline int numberL() const
234 {
235 return numberL_;
236 }
237
239 inline int baseL() const
240 {
241 return baseL_;
242 }
244 inline int maximumRowsExtra() const
245 {
246 return maximumRowsExtra_;
247 }
249 inline int numberColumns() const
250 {
251 return numberColumns_;
252 }
254 inline int numberElements() const
255 {
256 return totalElements_;
257 }
259 inline int numberForrestTomlin() const
260 {
262 }
264 inline int numberGoodColumns() const
265 {
266 return numberGoodU_;
267 }
269 inline double areaFactor() const
270 {
271 return areaFactor_;
272 }
273 inline void areaFactor(double value)
274 {
275 areaFactor_ = value;
276 }
278 double adjustedAreaFactor() const;
280 inline void relaxAccuracyCheck(double value)
281 {
282 relaxCheck_ = value;
283 }
284 inline double getAccuracyCheck() const
285 {
286 return relaxCheck_;
287 }
289 inline int messageLevel() const
290 {
291 return messageLevel_;
292 }
293 void messageLevel(int value);
295 inline int maximumPivots() const
296 {
297 return maximumPivots_;
298 }
299 void maximumPivots(int value);
300
302 inline int denseThreshold() const
303 {
304 return denseThreshold_;
305 }
307 inline void setDenseThreshold(int value)
308 {
309 denseThreshold_ = value;
310 }
312 inline double pivotTolerance() const
313 {
314 return pivotTolerance_;
315 }
316 void pivotTolerance(double value);
318 inline double zeroTolerance() const
319 {
320 return zeroTolerance_;
321 }
322 void zeroTolerance(double value);
323#ifndef COIN_FAST_CODE
325 inline double slackValue() const
326 {
327 return slackValue_;
328 }
329 void slackValue(double value);
330#endif
332 double maximumCoefficient() const;
334 inline bool forrestTomlin() const
335 {
336 return doForrestTomlin_;
337 }
338 inline void setForrestTomlin(bool value)
339 {
340 doForrestTomlin_ = value;
341 }
343 inline bool spaceForForrestTomlin() const
344 {
346 int space = lengthAreaU_ - (start + numberRowsExtra_);
347 return (space >= 0) && doForrestTomlin_;
348 }
350
353
355 inline int numberDense() const
356 {
357 return numberDense_;
358 }
359
361 inline int numberElementsU() const
362 {
363 return lengthU_;
364 }
366 inline void setNumberElementsU(int value)
367 {
368 lengthU_ = value;
369 }
371 inline int lengthAreaU() const
372 {
373 return lengthAreaU_;
374 }
376 inline int numberElementsL() const
377 {
378 return lengthL_;
379 }
381 inline int lengthAreaL() const
382 {
383 return lengthAreaL_;
384 }
386 inline int numberElementsR() const
387 {
388 return lengthR_;
389 }
391 inline int numberCompressions() const
392 {
393 return numberCompressions_;
394 }
396 inline int *numberInRow() const
397 {
398 return numberInRow_.array();
399 }
401 inline int *numberInColumn() const
402 {
403 return numberInColumn_.array();
404 }
407 {
408 return elementU_.array();
409 }
411 inline int *indexRowU() const
412 {
413 return indexRowU_.array();
414 }
416 inline int *startColumnU() const
417 {
418 return startColumnU_.array();
419 }
422 {
424 }
428 inline int biasLU() const
429 {
430 return biasLU_;
431 }
432 inline void setBiasLU(int value)
433 {
434 biasLU_ = value;
435 }
441 inline int persistenceFlag() const
442 {
443 return persistenceFlag_;
444 }
445 void setPersistenceFlag(int value);
447
450
459 int pivotRow,
460 double pivotCheck,
461 bool checkBeforeModifying = false,
462 double acceptablePivot = 1.0e-8);
468 int *deleted,
469 int internalPivotRow);
470#ifdef ABC_USE_COIN_FACTORIZATION
473 CoinIndexedVector *fakeVector(CoinIndexedVector *vector,
474 int already = 0) const;
475 void deleteFakeVector(CoinIndexedVector *vector,
476 CoinIndexedVector *fakeVector) const;
481 double checkReplacePart1(CoinIndexedVector *regionSparse,
482 int pivotRow);
487 double checkReplacePart1(CoinIndexedVector *regionSparse,
488 CoinIndexedVector *partialUpdate,
489 int pivotRow);
492 int checkReplacePart2(int pivotRow,
493 double btranAlpha,
494 double ftranAlpha,
495 double ftAlpha,
496 double acceptablePivot = 1.0e-8);
499 void replaceColumnPart3(CoinIndexedVector *regionSparse,
500 int pivotRow,
501 double alpha);
504 void replaceColumnPart3(CoinIndexedVector *regionSparse,
505 CoinIndexedVector *partialUpdate,
506 int pivotRow,
507 double alpha);
515 int updateColumnFT(CoinIndexedVector &regionSparse);
516 int updateColumnFTPart1(CoinIndexedVector &regionSparse);
517 void updateColumnFTPart2(CoinIndexedVector &regionSparse);
521 void updateColumnFT(CoinIndexedVector &regionSparseFT,
522 CoinIndexedVector &partialUpdate,
523 int which);
525 int updateColumn(CoinIndexedVector &regionSparse) const;
530 int updateTwoColumnsFT(CoinIndexedVector &regionSparseFT,
531 CoinIndexedVector &regionSparseOther);
533 int updateColumnTranspose(CoinIndexedVector &regionSparse) const;
535 void updateColumnCpu(CoinIndexedVector &regionSparse, int whichCpu) const;
537 void updateColumnTransposeCpu(CoinIndexedVector &regionSparse, int whichCpu) const;
539 void updateFullColumn(CoinIndexedVector &regionSparse) const;
541 void updateFullColumnTranspose(CoinIndexedVector &regionSparse) const;
543 void updateWeights(CoinIndexedVector &regionSparse) const;
545 inline bool wantsTableauColumn() const
546 {
547 return false;
548 }
550 inline double minimumPivotTolerance() const
551 {
552 return pivotTolerance_;
553 }
554 inline void minimumPivotTolerance(double value)
555 {
556 pivotTolerance(value);
557 }
559 inline void setParallelMode(int value)
560 {
561 parallelMode_ = value;
562 }
564 inline void setSolveMode(int value)
565 {
566 parallelMode_ &= 3;
567 parallelMode_ |= (value << 2);
568 }
570 inline int solveMode() const
571 {
572 return parallelMode_ >> 2;
573 }
575 void updatePartialUpdate(CoinIndexedVector &partialUpdate);
577 void makeNonSingular(int *COIN_RESTRICT sequence);
578#endif
579#if ABOCA_LITE_FACTORIZATION
581 void replaceColumn1(CoinIndexedVector *regionSparse, int pivotRow);
583 int replaceColumn2(CoinIndexedVector *regionSparse,
584 int pivotRow,
585 double pivotCheck);
586#endif
588
599 CoinIndexedVector *regionSparse2);
603 CoinIndexedVector *regionSparse2,
604 bool noPermute = false) const;
611 CoinIndexedVector *regionSparse2,
612 CoinIndexedVector *regionSparse3,
613 bool noPermuteRegion3 = false);
619 CoinIndexedVector *regionSparse2) const;
621 void updateOneColumnTranspose(CoinIndexedVector *regionWork, int &statistics) const;
627 CoinIndexedVector *regionSparse2,
628 CoinIndexedVector *regionSparse3,
629 int type) const;
631 void goSparse();
633 inline int sparseThreshold() const
634 {
635 return sparseThreshold_;
636 }
638 void sparseThreshold(int value);
640
641
646 inline void clearArrays()
647 {
649 }
651
654
658 int indicesRow[],
659 int indicesColumn[], double elements[]);
660
664 int indicesRow[], double elements[]);
665
669 int indicesColumn[], double elements[]);
670
672 int deleteColumn(int Row);
674 int deleteRow(int Row);
675
679 int replaceRow(int whichRow, int numberElements,
680 const int indicesColumn[], const double elements[]);
682 void emptyRows(int numberToEmpty, const int which[]);
684
688#if 0 //def CLP_FACTORIZATION_INSTRUMENT
689 inline bool collectStatistics() const
690 { return collectStatistics_;}
692 inline void setCollectStatistics(bool onOff) const
693 { collectStatistics_ = onOff;}
694#else
695 inline bool collectStatistics() const
696 {
697 return true;
698 }
700 inline void setCollectStatistics(bool onOff) const
701 {
702 }
703#endif
705 void gutsOfDestructor(int type = 1);
707 void gutsOfInitialize(int type);
708 void gutsOfCopy(const CoinFactorization &other);
709
712
714
718 int numberColumns,
719 int maximumL,
720 int maximumU);
721
724 void preProcess(int state,
725 int possibleDuplicates = -1);
727 int factor();
728
729protected:
742
744 bool pivotOneOtherRow(int pivotRow,
745 int pivotColumn);
747 bool pivotRowSingleton(int pivotRow,
748 int pivotColumn);
750 bool pivotColumnSingleton(int pivotRow,
751 int pivotColumn);
752
757 bool getColumnSpace(int iColumn,
758 int extraNeeded);
759
762 bool reorderU();
766 bool getColumnSpaceIterateR(int iColumn, double value,
767 int iRow);
773 int getColumnSpaceIterate(int iColumn, double value,
774 int iRow);
778 bool getRowSpace(int iRow, int extraNeeded);
779
783 bool getRowSpaceIterate(int iRow,
784 int extraNeeded);
788 inline void addLink(int index, int count)
789 {
790 int *nextCount = nextCount_.array();
791 int *firstCount = firstCount_.array();
792 int *lastCount = lastCount_.array();
793 int next = firstCount[count];
794 lastCount[index] = -2 - count;
795 if (next < 0) {
796 //first with that count
797 firstCount[count] = index;
798 nextCount[index] = -1;
799 } else {
800 firstCount[count] = index;
801 nextCount[index] = next;
802 lastCount[next] = index;
803 }
804 }
806 inline void deleteLink(int index)
807 {
808 int *nextCount = nextCount_.array();
809 int *firstCount = firstCount_.array();
810 int *lastCount = lastCount_.array();
811 int next = nextCount[index];
812 int last = lastCount[index];
813 if (last >= 0) {
814 nextCount[last] = next;
815 } else {
816 int count = -last - 2;
817
818 firstCount[count] = next;
819 }
820 if (next >= 0) {
821 lastCount[next] = last;
822 }
823 nextCount[index] = -2;
824 lastCount[index] = -2;
825 return;
826 }
828 void separateLinks(int count, bool rowsFirst);
830 void cleanup();
831
833 void updateColumnL(CoinIndexedVector *region, int *indexIn) const;
835 void updateColumnLDensish(CoinIndexedVector *region, int *indexIn) const;
837 void updateColumnLSparse(CoinIndexedVector *region, int *indexIn) const;
839 void updateColumnLSparsish(CoinIndexedVector *region, int *indexIn) const;
840
842 void updateColumnR(CoinIndexedVector *region) const;
845 void updateColumnRFT(CoinIndexedVector *region, int *indexIn);
846
848 void updateColumnU(CoinIndexedVector *region, int *indexIn) const;
849
852 int *indexIn) const;
855 int *indexIn) const;
858 int *COIN_RESTRICT regionIndex) const;
861 int &numberNonZero1,
862 double *COIN_RESTRICT region1,
863 int *COIN_RESTRICT index1,
864 int &numberNonZero2,
865 double *COIN_RESTRICT region2,
866 int *COIN_RESTRICT index2) const;
868 void updateColumnPFI(CoinIndexedVector *regionSparse) const;
870 void permuteBack(CoinIndexedVector *regionSparse,
871 CoinIndexedVector *outVector) const;
872
878 int smallestIndex) const;
882 int smallestIndex) const;
886 int smallestIndex) const;
893 int smallestIndex) const;
894
901
912
913public:
919 int pivotRow, double alpha);
920
921protected:
924 int checkPivot(double saveFromU, double oldPivot) const;
925 /********************************* START LARGE TEMPLATE ********/
926#ifdef INT_IS_8
927#define COINFACTORIZATION_BITS_PER_INT 64
928#define COINFACTORIZATION_SHIFT_PER_INT 6
929#define COINFACTORIZATION_MASK_PER_INT 0x3f
930#else
931#define COINFACTORIZATION_BITS_PER_INT 32
932#define COINFACTORIZATION_SHIFT_PER_INT 5
933#define COINFACTORIZATION_MASK_PER_INT 0x1f
934#endif
935 template < class T >
936 inline bool
937 pivot(int pivotRow,
938 int pivotColumn,
939 int pivotRowPosition,
940 int pivotColumnPosition,
942 unsigned int workArea2[],
943 int increment2,
944 T markRow[],
945 int largeInteger)
946 {
947 int *indexColumnU = indexColumnU_.array();
951 int *indexRowU = indexRowU_.array();
952 int *startRowU = startRowU_.array();
955 int *indexRowL = indexRowL_.array();
956 int *saveColumn = saveColumn_.array();
957 int *nextRow = nextRow_.array();
958 int *lastRow = lastRow_.array();
959
960 //store pivot columns (so can easily compress)
961 int numberInPivotRow = numberInRow[pivotRow] - 1;
962 int startColumn = startColumnU[pivotColumn];
963 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
964 int endColumn = startColumn + numberInPivotColumn + 1;
965 int put = 0;
966 int startRow = startRowU[pivotRow];
967 int endRow = startRow + numberInPivotRow + 1;
968
969 if (pivotColumnPosition < 0) {
970 for (pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++) {
971 int iColumn = indexColumnU[pivotColumnPosition];
972 if (iColumn != pivotColumn) {
973 saveColumn[put++] = iColumn;
974 } else {
975 break;
976 }
977 }
978 } else {
979 for (int i = startRow; i < pivotColumnPosition; i++) {
980 saveColumn[put++] = indexColumnU[i];
981 }
982 }
983 assert(pivotColumnPosition < endRow);
984 assert(indexColumnU[pivotColumnPosition] == pivotColumn);
985 pivotColumnPosition++;
986 for (; pivotColumnPosition < endRow; pivotColumnPosition++) {
987 saveColumn[put++] = indexColumnU[pivotColumnPosition];
988 }
989 //take out this bit of indexColumnU
990 int next = nextRow[pivotRow];
991 int last = lastRow[pivotRow];
992
993 nextRow[last] = next;
994 lastRow[next] = last;
995 nextRow[pivotRow] = numberGoodU_; //use for permute
996 lastRow[pivotRow] = -2;
997 numberInRow[pivotRow] = 0;
998 //store column in L, compress in U and take column out
999 int l = lengthL_;
1000
1001 if (l + numberInPivotColumn > lengthAreaL_) {
1002 //need more memory
1003 if ((messageLevel_ & 4) != 0)
1004 printf("more memory needed in middle of invert\n");
1005 return false;
1006 }
1007 //l+=currentAreaL_->elementByColumn-elementL;
1008 int lSave = l;
1009
1011 startColumnL[numberGoodL_] = l; //for luck and first time
1012 numberGoodL_++;
1013 startColumnL[numberGoodL_] = l + numberInPivotColumn;
1014 lengthL_ += numberInPivotColumn;
1015 if (pivotRowPosition < 0) {
1016 for (pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++) {
1017 int iRow = indexRowU[pivotRowPosition];
1018 if (iRow != pivotRow) {
1019 indexRowL[l] = iRow;
1020 elementL[l] = elementU[pivotRowPosition];
1021 markRow[iRow] = static_cast< T >(l - lSave);
1022 l++;
1023 //take out of row list
1024 int start = startRowU[iRow];
1025 int end = start + numberInRow[iRow];
1026 int where = start;
1027
1028 while (indexColumnU[where] != pivotColumn) {
1029 where++;
1030 } /* endwhile */
1031#if DEBUG_COIN
1032 if (where >= end) {
1033 abort();
1034 }
1035#endif
1036 indexColumnU[where] = indexColumnU[end - 1];
1037 numberInRow[iRow]--;
1038 } else {
1039 break;
1040 }
1041 }
1042 } else {
1043 int i;
1044
1045 for (i = startColumn; i < pivotRowPosition; i++) {
1046 int iRow = indexRowU[i];
1047
1048 markRow[iRow] = static_cast< T >(l - lSave);
1049 indexRowL[l] = iRow;
1050 elementL[l] = elementU[i];
1051 l++;
1052 //take out of row list
1053 int start = startRowU[iRow];
1054 int end = start + numberInRow[iRow];
1055 int where = start;
1056
1057 while (indexColumnU[where] != pivotColumn) {
1058 where++;
1059 } /* endwhile */
1060#if DEBUG_COIN
1061 if (where >= end) {
1062 abort();
1063 }
1064#endif
1065 indexColumnU[where] = indexColumnU[end - 1];
1066 numberInRow[iRow]--;
1067 assert(numberInRow[iRow] >= 0);
1068 }
1069 }
1070 assert(pivotRowPosition < endColumn);
1071 assert(indexRowU[pivotRowPosition] == pivotRow);
1072 CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
1073 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1074
1075 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1076 pivotRowPosition++;
1077 for (; pivotRowPosition < endColumn; pivotRowPosition++) {
1078 int iRow = indexRowU[pivotRowPosition];
1079
1080 markRow[iRow] = static_cast< T >(l - lSave);
1081 indexRowL[l] = iRow;
1082 elementL[l] = elementU[pivotRowPosition];
1083 l++;
1084 //take out of row list
1085 int start = startRowU[iRow];
1086 int end = start + numberInRow[iRow];
1087 int where = start;
1088
1089 while (indexColumnU[where] != pivotColumn) {
1090 where++;
1091 } /* endwhile */
1092#if DEBUG_COIN
1093 if (where >= end) {
1094 abort();
1095 }
1096#endif
1097 indexColumnU[where] = indexColumnU[end - 1];
1098 numberInRow[iRow]--;
1099 assert(numberInRow[iRow] >= 0);
1100 }
1101 markRow[pivotRow] = static_cast< T >(largeInteger);
1102 //compress pivot column (move pivot to front including saved)
1104 //use end of L for temporary space
1105 int *indexL = &indexRowL[lSave];
1106 CoinFactorizationDouble *multipliersL = &elementL[lSave];
1107
1108 //adjust
1109 int j;
1110
1111 for (j = 0; j < numberInPivotColumn; j++) {
1112 multipliersL[j] *= pivotMultiplier;
1113 }
1114 //zero out fill
1115 int iErase;
1116 for (iErase = 0; iErase < increment2 * numberInPivotRow;
1117 iErase++) {
1118 workArea2[iErase] = 0;
1119 }
1120 int added = numberInPivotRow * numberInPivotColumn;
1121 unsigned int *temp2 = workArea2;
1122 int *nextColumn = nextColumn_.array();
1123
1124 //pack down and move to work
1125 int jColumn;
1126 for (jColumn = 0; jColumn < numberInPivotRow; jColumn++) {
1127 int iColumn = saveColumn[jColumn];
1128 int startColumn = startColumnU[iColumn];
1129 int endColumn = startColumn + numberInColumn[iColumn];
1130 int iRow = indexRowU[startColumn];
1131 CoinFactorizationDouble value = elementU[startColumn];
1132 double largest;
1133 int put = startColumn;
1134 int positionLargest = -1;
1135 CoinFactorizationDouble thisPivotValue = 0.0;
1136
1137 //compress column and find largest not updated
1138 bool checkLargest;
1139 int mark = markRow[iRow];
1140
1141 if (mark == largeInteger + 1) {
1142 largest = fabs(value);
1143 positionLargest = put;
1144 put++;
1145 checkLargest = false;
1146 } else {
1147 //need to find largest
1148 largest = 0.0;
1149 checkLargest = true;
1150 if (mark != largeInteger) {
1151 //will be updated
1152 work[mark] = value;
1153 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1154 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1155
1156 temp2[word] = temp2[word] | (1 << bit); //say already in counts
1157 added--;
1158 } else {
1159 thisPivotValue = value;
1160 }
1161 }
1162 int i;
1163 for (i = startColumn + 1; i < endColumn; i++) {
1164 iRow = indexRowU[i];
1165 value = elementU[i];
1166 int mark = markRow[iRow];
1167
1168 if (mark == largeInteger + 1) {
1169 //keep
1170 indexRowU[put] = iRow;
1171 elementU[put] = value;
1172 if (checkLargest) {
1173 double absValue = fabs(value);
1174
1175 if (absValue > largest) {
1176 largest = absValue;
1177 positionLargest = put;
1178 }
1179 }
1180 put++;
1181 } else if (mark != largeInteger) {
1182 //will be updated
1183 work[mark] = value;
1184 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1185 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1186
1187 temp2[word] = temp2[word] | (1 << bit); //say already in counts
1188 added--;
1189 } else {
1190 thisPivotValue = value;
1191 }
1192 }
1193 //slot in pivot
1194 elementU[put] = elementU[startColumn];
1195 indexRowU[put] = indexRowU[startColumn];
1196 if (positionLargest == startColumn) {
1197 positionLargest = put; //follow if was largest
1198 }
1199 put++;
1200 elementU[startColumn] = thisPivotValue;
1201 indexRowU[startColumn] = pivotRow;
1202 //clean up counts
1203 startColumn++;
1204 numberInColumn[iColumn] = put - startColumn;
1205 int *numberInColumnPlus = numberInColumnPlus_.array();
1206 numberInColumnPlus[iColumn]++;
1207 startColumnU[iColumn]++;
1208 //how much space have we got
1209 int next = nextColumn[iColumn];
1210 int space;
1211
1212 space = startColumnU[next] - put - numberInColumnPlus[next];
1213 //assume no zero elements
1214 if (numberInPivotColumn > space) {
1215 //getColumnSpace also moves fixed part
1216 if (!getColumnSpace(iColumn, numberInPivotColumn)) {
1217 return false;
1218 }
1219 //redo starts
1220 if (positionLargest >= 0)
1221 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
1222 startColumn = startColumnU[iColumn];
1223 put = startColumn + numberInColumn[iColumn];
1224 }
1225 double tolerance = zeroTolerance_;
1226
1227 int *nextCount = nextCount_.array();
1228 for (j = 0; j < numberInPivotColumn; j++) {
1229 value = work[j] - thisPivotValue * multipliersL[j];
1230 double absValue = fabs(value);
1231
1232 if (absValue > tolerance) {
1233 work[j] = 0.0;
1234 assert(put < lengthAreaU_);
1235 elementU[put] = value;
1236 indexRowU[put] = indexL[j];
1237 if (absValue > largest) {
1238 largest = absValue;
1239 positionLargest = put;
1240 }
1241 put++;
1242 } else {
1243 work[j] = 0.0;
1244 added--;
1245 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1246 int bit = j & COINFACTORIZATION_MASK_PER_INT;
1247
1248 if (temp2[word] & (1 << bit)) {
1249 //take out of row list
1250 iRow = indexL[j];
1251 int start = startRowU[iRow];
1252 int end = start + numberInRow[iRow];
1253 int where = start;
1254
1255 while (indexColumnU[where] != iColumn) {
1256 where++;
1257 } /* endwhile */
1258#if DEBUG_COIN
1259 if (where >= end) {
1260 abort();
1261 }
1262#endif
1263 indexColumnU[where] = indexColumnU[end - 1];
1264 numberInRow[iRow]--;
1265 } else {
1266 //make sure won't be added
1267 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1268 int bit = j & COINFACTORIZATION_MASK_PER_INT;
1269
1270 temp2[word] = temp2[word] | (1 << bit); //say already in counts
1271 }
1272 }
1273 }
1274 numberInColumn[iColumn] = put - startColumn;
1275 //move largest
1276 if (positionLargest >= 0) {
1277 value = elementU[positionLargest];
1278 iRow = indexRowU[positionLargest];
1279 elementU[positionLargest] = elementU[startColumn];
1280 indexRowU[positionLargest] = indexRowU[startColumn];
1281 elementU[startColumn] = value;
1282 indexRowU[startColumn] = iRow;
1283 }
1284 //linked list for column
1285 if (nextCount[iColumn + numberRows_] != -2) {
1286 //modify linked list
1287 deleteLink(iColumn + numberRows_);
1288 addLink(iColumn + numberRows_, numberInColumn[iColumn]);
1289 }
1290 temp2 += increment2;
1291 }
1292 //get space for row list
1293 unsigned int *putBase = workArea2;
1294 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1295 int i = 0;
1296
1297 // do linked lists and update counts
1298 while (bigLoops) {
1299 bigLoops--;
1300 int bit;
1301 for (bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++) {
1302 unsigned int *putThis = putBase;
1303 int iRow = indexL[i];
1304
1305 //get space
1306 int number = 0;
1307 int jColumn;
1308
1309 for (jColumn = 0; jColumn < numberInPivotRow; jColumn++) {
1310 unsigned int test = *putThis;
1311
1312 putThis += increment2;
1313 test = 1 - ((test >> bit) & 1);
1314 number += test;
1315 }
1316 int next = nextRow[iRow];
1317 int space;
1318
1319 space = startRowU[next] - startRowU[iRow];
1320 number += numberInRow[iRow];
1321 if (space < number) {
1322 if (!getRowSpace(iRow, number)) {
1323 return false;
1324 }
1325 }
1326 // now do
1327 putThis = putBase;
1328 next = nextRow[iRow];
1329 number = numberInRow[iRow];
1330 int end = startRowU[iRow] + number;
1331 int saveIndex = indexColumnU[startRowU[next]];
1332
1333 //add in
1334 for (jColumn = 0; jColumn < numberInPivotRow; jColumn++) {
1335 unsigned int test = *putThis;
1336
1337 putThis += increment2;
1338 test = 1 - ((test >> bit) & 1);
1339 indexColumnU[end] = saveColumn[jColumn];
1340 end += test;
1341 }
1342 //put back next one in case zapped
1343 indexColumnU[startRowU[next]] = saveIndex;
1344 markRow[iRow] = static_cast< T >(largeInteger + 1);
1345 number = end - startRowU[iRow];
1346 numberInRow[iRow] = number;
1347 deleteLink(iRow);
1348 addLink(iRow, number);
1349 }
1350 putBase++;
1351 } /* endwhile */
1352 int bit;
1353
1354 for (bit = 0; i < numberInPivotColumn; i++, bit++) {
1355 unsigned int *putThis = putBase;
1356 int iRow = indexL[i];
1357
1358 //get space
1359 int number = 0;
1360 int jColumn;
1361
1362 for (jColumn = 0; jColumn < numberInPivotRow; jColumn++) {
1363 unsigned int test = *putThis;
1364
1365 putThis += increment2;
1366 test = 1 - ((test >> bit) & 1);
1367 number += test;
1368 }
1369 int next = nextRow[iRow];
1370 int space;
1371
1372 space = startRowU[next] - startRowU[iRow];
1373 number += numberInRow[iRow];
1374 if (space < number) {
1375 if (!getRowSpace(iRow, number)) {
1376 return false;
1377 }
1378 }
1379 // now do
1380 putThis = putBase;
1381 next = nextRow[iRow];
1382 number = numberInRow[iRow];
1383 int end = startRowU[iRow] + number;
1384 int saveIndex;
1385
1386 saveIndex = indexColumnU[startRowU[next]];
1387
1388 //add in
1389 for (jColumn = 0; jColumn < numberInPivotRow; jColumn++) {
1390 unsigned int test = *putThis;
1391
1392 putThis += increment2;
1393 test = 1 - ((test >> bit) & 1);
1394
1395 indexColumnU[end] = saveColumn[jColumn];
1396 end += test;
1397 }
1398 indexColumnU[startRowU[next]] = saveIndex;
1399 markRow[iRow] = static_cast< T >(largeInteger + 1);
1400 number = end - startRowU[iRow];
1401 numberInRow[iRow] = number;
1402 deleteLink(iRow);
1403 addLink(iRow, number);
1404 }
1405 markRow[pivotRow] = static_cast< T >(largeInteger + 1);
1406 //modify linked list for pivots
1407 deleteLink(pivotRow);
1409 totalElements_ += added;
1410 return true;
1411 }
1412
1413 /********************************* END LARGE TEMPLATE ********/
1415
1416protected:
1423#ifndef COIN_FAST_CODE
1426#else
1427#ifndef slackValue_
1428#define slackValue_ -1.0
1429#endif
1430#endif
1470
1475 //int increasingRows_;
1476
1481
1484
1487
1490
1494
1497
1500
1503
1506
1509
1512
1515
1518
1521
1524
1527
1530
1533
1536
1539
1542
1544 //int baseU_;
1545
1548
1551
1554
1557
1560
1563
1566
1569
1572
1575
1578
1581
1584
1587
1590
1593
1596
1599
1602
1605
1607 double *denseArea_;
1608
1611
1614
1617
1620
1623
1626
1629
1630public:
1632 mutable double ftranCountInput_;
1633 mutable double ftranCountAfterL_;
1634 mutable double ftranCountAfterR_;
1635 mutable double ftranCountAfterU_;
1636 mutable double btranCountInput_;
1637 mutable double btranCountAfterU_;
1638 mutable double btranCountAfterR_;
1639 mutable double btranCountAfterL_;
1640
1644
1652
1653protected:
1655#if 0
1656 mutable bool collectStatistics_;
1657#else
1658#define collectStatistics_ 1
1659#endif
1660
1663
1666
1669
1672
1675
1678#if ABOCA_LITE_FACTORIZATION
1680 int sparseOffset_;
1681#endif
1692#ifdef ABC_USE_COIN_FACTORIZATION
1694 int parallelMode_;
1695#endif
1697};
1698// Dense coding
1699#ifdef INTEL_COMPILER
1700#define COIN_FACTORIZATION_DENSE_CODE 3
1701#endif
1702#ifdef COIN_HAS_LAPACK
1703#ifndef COIN_FACTORIZATION_DENSE_CODE
1704#define COIN_FACTORIZATION_DENSE_CODE 1
1705#endif
1706#endif
1707#ifdef COIN_FACTORIZATION_DENSE_CODE
1708/* Type of Fortran integer translated into C */
1709#ifndef ipfint
1710//typedef ipfint FORTRAN_INTEGER_TYPE ;
1711typedef int ipfint;
1712typedef const int cipfint;
1713#endif
1714#endif
1715#endif
1716// Extra for ugly include
1717#ifdef UGLY_COIN_FACTOR_CODING
1718#define FAC_UNSET (FAC_SET + 1)
1719{
1720 goodPivot = false;
1721 //store pivot columns (so can easily compress)
1722 int startColumnThis = startColumn[iPivotColumn];
1723 int endColumn = startColumnThis + numberDoColumn + 1;
1724 int put = 0;
1725 int startRowThis = startRow[iPivotRow];
1726 int endRow = startRowThis + numberDoRow + 1;
1727 if (pivotColumnPosition < 0) {
1728 for (pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++) {
1729 int iColumn = indexColumn[pivotColumnPosition];
1730 if (iColumn != iPivotColumn) {
1731 saveColumn[put++] = iColumn;
1732 } else {
1733 break;
1734 }
1735 }
1736 } else {
1737 for (int i = startRowThis; i < pivotColumnPosition; i++) {
1738 saveColumn[put++] = indexColumn[i];
1739 }
1740 }
1741 assert(pivotColumnPosition < endRow);
1742 assert(indexColumn[pivotColumnPosition] == iPivotColumn);
1743 pivotColumnPosition++;
1744 for (; pivotColumnPosition < endRow; pivotColumnPosition++) {
1745 saveColumn[put++] = indexColumn[pivotColumnPosition];
1746 }
1747 //take out this bit of indexColumn
1748 int next = nextRow[iPivotRow];
1749 int last = lastRow[iPivotRow];
1750
1751 nextRow[last] = next;
1752 lastRow[next] = last;
1753 nextRow[iPivotRow] = numberGoodU_; //use for permute
1754 lastRow[iPivotRow] = -2;
1755 numberInRow[iPivotRow] = 0;
1756 //store column in L, compress in U and take column out
1757 int l = lengthL_;
1758 // **** HORRID coding coming up but a goto seems best!
1759 {
1760 if (l + numberDoColumn > lengthAreaL_) {
1761 //need more memory
1762 if ((messageLevel_ & 4) != 0)
1763 printf("more memory needed in middle of invert\n");
1764 goto BAD_PIVOT;
1765 }
1766 //l+=currentAreaL_->elementByColumn-elementL;
1767 int lSave = l;
1768
1770 startColumnL[numberGoodL_] = l; //for luck and first time
1771 numberGoodL_++;
1772 startColumnL[numberGoodL_] = l + numberDoColumn;
1773 lengthL_ += numberDoColumn;
1774 if (pivotRowPosition < 0) {
1775 for (pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++) {
1776 int iRow = indexRow[pivotRowPosition];
1777 if (iRow != iPivotRow) {
1778 indexRowL[l] = iRow;
1779 elementL[l] = element[pivotRowPosition];
1780 markRow[iRow] = l - lSave;
1781 l++;
1782 //take out of row list
1783 int start = startRow[iRow];
1784 int end = start + numberInRow[iRow];
1785 int where = start;
1786
1787 while (indexColumn[where] != iPivotColumn) {
1788 where++;
1789 } /* endwhile */
1790#if DEBUG_COIN
1791 if (where >= end) {
1792 abort();
1793 }
1794#endif
1795 indexColumn[where] = indexColumn[end - 1];
1796 numberInRow[iRow]--;
1797 } else {
1798 break;
1799 }
1800 }
1801 } else {
1802 int i;
1803
1804 for (i = startColumnThis; i < pivotRowPosition; i++) {
1805 int iRow = indexRow[i];
1806
1807 markRow[iRow] = l - lSave;
1808 indexRowL[l] = iRow;
1809 elementL[l] = element[i];
1810 l++;
1811 //take out of row list
1812 int start = startRow[iRow];
1813 int end = start + numberInRow[iRow];
1814 int where = start;
1815
1816 while (indexColumn[where] != iPivotColumn) {
1817 where++;
1818 } /* endwhile */
1819#if DEBUG_COIN
1820 if (where >= end) {
1821 abort();
1822 }
1823#endif
1824 indexColumn[where] = indexColumn[end - 1];
1825 numberInRow[iRow]--;
1826 assert(numberInRow[iRow] >= 0);
1827 }
1828 }
1829 assert(pivotRowPosition < endColumn);
1830 assert(indexRow[pivotRowPosition] == iPivotRow);
1831 CoinFactorizationDouble pivotElement = element[pivotRowPosition];
1832 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1833
1834 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1835 pivotRowPosition++;
1836 for (; pivotRowPosition < endColumn; pivotRowPosition++) {
1837 int iRow = indexRow[pivotRowPosition];
1838
1839 markRow[iRow] = l - lSave;
1840 indexRowL[l] = iRow;
1841 elementL[l] = element[pivotRowPosition];
1842 l++;
1843 //take out of row list
1844 int start = startRow[iRow];
1845 int end = start + numberInRow[iRow];
1846 int where = start;
1847
1848 while (indexColumn[where] != iPivotColumn) {
1849 where++;
1850 } /* endwhile */
1851#if DEBUG_COIN
1852 if (where >= end) {
1853 abort();
1854 }
1855#endif
1856 indexColumn[where] = indexColumn[end - 1];
1857 numberInRow[iRow]--;
1858 assert(numberInRow[iRow] >= 0);
1859 }
1860 markRow[iPivotRow] = FAC_SET;
1861 //compress pivot column (move pivot to front including saved)
1862 numberInColumn[iPivotColumn] = 0;
1863 //use end of L for temporary space
1864 int *indexL = &indexRowL[lSave];
1865 CoinFactorizationDouble *multipliersL = &elementL[lSave];
1866
1867 //adjust
1868 int j;
1869
1870 for (j = 0; j < numberDoColumn; j++) {
1871 multipliersL[j] *= pivotMultiplier;
1872 }
1873 //zero out fill
1874 int iErase;
1875 for (iErase = 0; iErase < increment2 * numberDoRow;
1876 iErase++) {
1877 workArea2[iErase] = 0;
1878 }
1879 int added = numberDoRow * numberDoColumn;
1880 unsigned int *temp2 = workArea2;
1881 int *nextColumn = nextColumn_.array();
1882
1883 //pack down and move to work
1884 int jColumn;
1885 for (jColumn = 0; jColumn < numberDoRow; jColumn++) {
1886 int iColumn = saveColumn[jColumn];
1887 int startColumnThis = startColumn[iColumn];
1888 int endColumn = startColumnThis + numberInColumn[iColumn];
1889 int iRow = indexRow[startColumnThis];
1890 CoinFactorizationDouble value = element[startColumnThis];
1891 double largest;
1892 int put = startColumnThis;
1893 int positionLargest = -1;
1894 CoinFactorizationDouble thisPivotValue = 0.0;
1895
1896 //compress column and find largest not updated
1897 bool checkLargest;
1898 int mark = markRow[iRow];
1899
1900 if (mark == FAC_UNSET) {
1901 largest = fabs(value);
1902 positionLargest = put;
1903 put++;
1904 checkLargest = false;
1905 } else {
1906 //need to find largest
1907 largest = 0.0;
1908 checkLargest = true;
1909 if (mark != FAC_SET) {
1910 //will be updated
1911 workArea[mark] = value;
1912 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1913 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1914
1915 temp2[word] = temp2[word] | (1 << bit); //say already in counts
1916 added--;
1917 } else {
1918 thisPivotValue = value;
1919 }
1920 }
1921 int i;
1922 for (i = startColumnThis + 1; i < endColumn; i++) {
1923 iRow = indexRow[i];
1924 value = element[i];
1925 int mark = markRow[iRow];
1926
1927 if (mark == FAC_UNSET) {
1928 //keep
1929 indexRow[put] = iRow;
1930 element[put] = value;
1931 if (checkLargest) {
1932 double absValue = fabs(value);
1933
1934 if (absValue > largest) {
1935 largest = absValue;
1936 positionLargest = put;
1937 }
1938 }
1939 put++;
1940 } else if (mark != FAC_SET) {
1941 //will be updated
1942 workArea[mark] = value;
1943 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1944 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1945
1946 temp2[word] = temp2[word] | (1 << bit); //say already in counts
1947 added--;
1948 } else {
1949 thisPivotValue = value;
1950 }
1951 }
1952 //slot in pivot
1953 element[put] = element[startColumnThis];
1954 indexRow[put] = indexRow[startColumnThis];
1955 if (positionLargest == startColumnThis) {
1956 positionLargest = put; //follow if was largest
1957 }
1958 put++;
1959 element[startColumnThis] = thisPivotValue;
1960 indexRow[startColumnThis] = iPivotRow;
1961 //clean up counts
1962 startColumnThis++;
1963 numberInColumn[iColumn] = put - startColumnThis;
1964 int *numberInColumnPlus = numberInColumnPlus_.array();
1965 numberInColumnPlus[iColumn]++;
1966 startColumn[iColumn]++;
1967 //how much space have we got
1968 int next = nextColumn[iColumn];
1969 int space;
1970
1971 space = startColumn[next] - put - numberInColumnPlus[next];
1972 //assume no zero elements
1973 if (numberDoColumn > space) {
1974 //getColumnSpace also moves fixed part
1975 if (!getColumnSpace(iColumn, numberDoColumn)) {
1976 goto BAD_PIVOT;
1977 }
1978 //redo starts
1979 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
1980 startColumnThis = startColumn[iColumn];
1981 put = startColumnThis + numberInColumn[iColumn];
1982 }
1983 double tolerance = zeroTolerance_;
1984
1985 int *nextCount = nextCount_.array();
1986 for (j = 0; j < numberDoColumn; j++) {
1987 value = workArea[j] - thisPivotValue * multipliersL[j];
1988 double absValue = fabs(value);
1989
1990 if (absValue > tolerance) {
1991 workArea[j] = 0.0;
1992 element[put] = value;
1993 indexRow[put] = indexL[j];
1994 if (absValue > largest) {
1995 largest = absValue;
1996 positionLargest = put;
1997 }
1998 put++;
1999 } else {
2000 workArea[j] = 0.0;
2001 added--;
2002 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
2003 int bit = j & COINFACTORIZATION_MASK_PER_INT;
2004
2005 if (temp2[word] & (1 << bit)) {
2006 //take out of row list
2007 iRow = indexL[j];
2008 int start = startRow[iRow];
2009 int end = start + numberInRow[iRow];
2010 int where = start;
2011
2012 while (indexColumn[where] != iColumn) {
2013 where++;
2014 } /* endwhile */
2015#if DEBUG_COIN
2016 if (where >= end) {
2017 abort();
2018 }
2019#endif
2020 indexColumn[where] = indexColumn[end - 1];
2021 numberInRow[iRow]--;
2022 } else {
2023 //make sure won't be added
2024 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
2025 int bit = j & COINFACTORIZATION_MASK_PER_INT;
2026
2027 temp2[word] = temp2[word] | (1 << bit); //say already in counts
2028 }
2029 }
2030 }
2031 numberInColumn[iColumn] = put - startColumnThis;
2032 //move largest
2033 if (positionLargest >= 0) {
2034 value = element[positionLargest];
2035 iRow = indexRow[positionLargest];
2036 element[positionLargest] = element[startColumnThis];
2037 indexRow[positionLargest] = indexRow[startColumnThis];
2038 element[startColumnThis] = value;
2039 indexRow[startColumnThis] = iRow;
2040 }
2041 //linked list for column
2042 if (nextCount[iColumn + numberRows_] != -2) {
2043 //modify linked list
2044 deleteLink(iColumn + numberRows_);
2045 addLink(iColumn + numberRows_, numberInColumn[iColumn]);
2046 }
2047 temp2 += increment2;
2048 }
2049 //get space for row list
2050 unsigned int *putBase = workArea2;
2051 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
2052 int i = 0;
2053
2054 // do linked lists and update counts
2055 while (bigLoops) {
2056 bigLoops--;
2057 int bit;
2058 for (bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++) {
2059 unsigned int *putThis = putBase;
2060 int iRow = indexL[i];
2061
2062 //get space
2063 int number = 0;
2064 int jColumn;
2065
2066 for (jColumn = 0; jColumn < numberDoRow; jColumn++) {
2067 unsigned int test = *putThis;
2068
2069 putThis += increment2;
2070 test = 1 - ((test >> bit) & 1);
2071 number += test;
2072 }
2073 int next = nextRow[iRow];
2074 int space;
2075
2076 space = startRow[next] - startRow[iRow];
2077 number += numberInRow[iRow];
2078 if (space < number) {
2079 if (!getRowSpace(iRow, number)) {
2080 goto BAD_PIVOT;
2081 }
2082 }
2083 // now do
2084 putThis = putBase;
2085 next = nextRow[iRow];
2086 number = numberInRow[iRow];
2087 int end = startRow[iRow] + number;
2088 int saveIndex = indexColumn[startRow[next]];
2089
2090 //add in
2091 for (jColumn = 0; jColumn < numberDoRow; jColumn++) {
2092 unsigned int test = *putThis;
2093
2094 putThis += increment2;
2095 test = 1 - ((test >> bit) & 1);
2096 indexColumn[end] = saveColumn[jColumn];
2097 end += test;
2098 }
2099 //put back next one in case zapped
2100 indexColumn[startRow[next]] = saveIndex;
2101 markRow[iRow] = FAC_UNSET;
2102 number = end - startRow[iRow];
2103 numberInRow[iRow] = number;
2104 deleteLink(iRow);
2105 addLink(iRow, number);
2106 }
2107 putBase++;
2108 } /* endwhile */
2109 int bit;
2110
2111 for (bit = 0; i < numberDoColumn; i++, bit++) {
2112 unsigned int *putThis = putBase;
2113 int iRow = indexL[i];
2114
2115 //get space
2116 int number = 0;
2117 int jColumn;
2118
2119 for (jColumn = 0; jColumn < numberDoRow; jColumn++) {
2120 unsigned int test = *putThis;
2121
2122 putThis += increment2;
2123 test = 1 - ((test >> bit) & 1);
2124 number += test;
2125 }
2126 int next = nextRow[iRow];
2127 int space;
2128
2129 space = startRow[next] - startRow[iRow];
2130 number += numberInRow[iRow];
2131 if (space < number) {
2132 if (!getRowSpace(iRow, number)) {
2133 goto BAD_PIVOT;
2134 }
2135 }
2136 // now do
2137 putThis = putBase;
2138 next = nextRow[iRow];
2139 number = numberInRow[iRow];
2140 int end = startRow[iRow] + number;
2141 int saveIndex;
2142
2143 saveIndex = indexColumn[startRow[next]];
2144
2145 //add in
2146 for (jColumn = 0; jColumn < numberDoRow; jColumn++) {
2147 unsigned int test = *putThis;
2148
2149 putThis += increment2;
2150 test = 1 - ((test >> bit) & 1);
2151
2152 indexColumn[end] = saveColumn[jColumn];
2153 end += test;
2154 }
2155 indexColumn[startRow[next]] = saveIndex;
2156 markRow[iRow] = FAC_UNSET;
2157 number = end - startRow[iRow];
2158 numberInRow[iRow] = number;
2159 deleteLink(iRow);
2160 addLink(iRow, number);
2161 }
2162 markRow[iPivotRow] = FAC_UNSET;
2163 //modify linked list for pivots
2164 deleteLink(iPivotRow);
2165 deleteLink(iPivotColumn + numberRows_);
2166 totalElements_ += added;
2167 goodPivot = true;
2168 // **** UGLY UGLY UGLY
2169 }
2170BAD_PIVOT:
2171
2172 ;
2173}
2174#undef FAC_UNSET
2175#endif
2176
2177/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2178*/
#define collectStatistics_
For statistics.
#define COINFACTORIZATION_BITS_PER_INT
#define COINFACTORIZATION_MASK_PER_INT
#define COINFACTORIZATION_SHIFT_PER_INT
#define COIN_RESTRICT
double CoinFactorizationDouble
Definition: CoinTypes.hpp:57
CoinFactorizationDouble * version.
CoinFactorizationDouble * array() const
Get Array.
This deals with Factorization and Updates.
int replaceColumn(CoinIndexedVector *regionSparse, int pivotRow, double pivotCheck, bool checkBeforeModifying=false, double acceptablePivot=1.0e-8)
Replaces one Column to basis, returns 0=OK, 1=Probably OK, 2=singular, 3=no room If checkBeforeModify...
int lengthAreaL_
Length of area reserved for L.
void updateColumnL(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL)
void getAreas(int numberRows, int numberColumns, int maximumL, int maximumU)
Gets space for a factorization, called by constructors.
void updateColumnLSparse(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when sparse.
double ftranAverageAfterL_
While these are average ratios collected over last period.
int numberGoodL_
Number factorized in L.
int lengthAreaU_
Length of area reserved for U.
int numberElementsL() const
Returns number in L area.
CoinFactorization()
Default constructor.
int baseL() const
Base of L.
int numberRows() const
Number of Rows after factorization.
void cleanup()
Cleans up at end of factorization.
CoinIntArrayWithLength pivotRowL_
Pivots for L.
int restoreFactorization(const char *file, bool factor=false)
Debug - restore from file - 0 if no error on file.
bool forrestTomlin() const
true if Forrest Tomlin update, false if PFI
void updateColumnTransposeUDensish(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) when densish, assumes index is sorted i....
void updateColumnTransposeLSparsish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparsish by row.
int maximumColumnsExtra_
Maximum number of Columns after iterating.
int messageLevel_
Detail in messages.
CoinIntArrayWithLength startColumnL_
Start of each column in L.
void setNumberElementsU(int value)
Setss number in U area.
CoinUnsignedIntArrayWithLength workArea2_
Second work area.
void maximumPivots(int value)
CoinIntArrayWithLength indexColumnL_
Index of column in row for L.
void updateColumnLDensish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when densish.
~CoinFactorization()
Destructor.
CoinFactorizationDoubleArrayWithLength elementU_
Elements of U.
CoinIntArrayWithLength sparse_
Sparse regions.
void pivotTolerance(double value)
int maximumPivots_
Maximum number of pivots before factorization.
double adjustedAreaFactor() const
Returns areaFactor but adjusted for dense.
int numberDense_
Number of dense rows.
int status() const
Returns status.
void updateColumnTransposeU(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU), assumes index is sorted i.e.
void messageLevel(int value)
void zeroTolerance(double value)
int numberRowsExtra_
Number of Rows after iterating.
void updateColumnTransposeR(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR)
int numberElements() const
Total number of elements in factorization.
void permuteBack(CoinIndexedVector *regionSparse, CoinIndexedVector *outVector) const
Permutes back at end of updateColumn.
void updateColumnTransposeUSparsish(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) when sparsish, assumes index is sorted i....
int maximumColumnsExtra()
Maximum number of Columns after iterating.
int numberL_
Number in L.
double areaFactor_
How much to multiply areas by.
void updateColumnTransposeUSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANU) when sparse, assumes index is sorted i.e.
CoinIntArrayWithLength startRowL_
Start of each row in L.
double * denseArea_
Dense area.
void show_self() const
Debug show object (shows one representation)
int * indexRowU() const
Row indices of U.
bool spaceForForrestTomlin() const
True if FT update and space.
CoinIntArrayWithLength saveColumn_
Columns left to do in a single pivot.
CoinIntArrayWithLength numberInColumn_
Number in each Column.
int getColumnSpaceIterate(int iColumn, double value, int iRow)
getColumnSpaceIterate.
void updateColumnTransposeLByRow(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when densish by row.
friend void CoinFactorizationUnitTest(const std::string &mpsDir)
void updateTwoColumnsUDensish(int &numberNonZero1, double *COIN_RESTRICT region1, int *COIN_RESTRICT index1, int &numberNonZero2, double *COIN_RESTRICT region2, int *COIN_RESTRICT index2) const
Updates part of 2 columns (FTRANU) real work.
int updateColumnTranspose(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2) const
Updates one column (BTRAN) from regionSparse2 regionSparse starts as zero and is zero at end Note - i...
CoinFactorizationDouble * elementR_
Elements of R.
void updateColumnU(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANU)
int numberGoodU_
Number factorized in U (not row singletons)
int checkPivot(double saveFromU, double oldPivot) const
Returns accuracy status of replaceColumn returns 0=OK, 1=Probably OK, 2=singular.
bool getRowSpaceIterate(int iRow, int extraNeeded)
Gets space for one Row with given length while iterating, may have to do compression (returns True if...
int * permuteBack() const
Returns address of permuteBack region.
int * pivotColumn() const
Returns address of pivotColumn region (also used for permuting)
void checkSparse()
See if worth going sparse.
double getAccuracyCheck() const
bool reorderU()
Reorders U so contiguous and in order (if there is space) Returns true if it could.
CoinIntArrayWithLength pivotColumn_
Pivot order for each Column.
CoinFactorizationDouble * elementByRowL() const
Elements in L (row copy)
void separateLinks(int count, bool rowsFirst)
Separate out links with same row/column count.
int maximumRowsExtra() const
Maximum of Rows after iterating.
int numberColumnsExtra_
Number of Columns after iterating.
CoinFactorizationDoubleArrayWithLength elementByRowL_
Elements in L (row copy)
int numberCompressions_
Number of compressions done.
void relaxAccuracyCheck(double value)
Allows change of pivot accuracy check 1.0 == none >1.0 relaxed.
CoinIntArrayWithLength pivotColumnBack_
Inverse Pivot order for each Column.
double ftranCountInput_
Below are all to collect.
void updateColumnUSparsish(CoinIndexedVector *regionSparse, int *indexIn) const
Updates part of column (FTRANU) when sparsish.
void updateColumnTransposeLSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparse (by Row)
CoinIntArrayWithLength markRow_
Marks rows to be updated.
int updateTwoColumnsFT(CoinIndexedVector *regionSparse1, CoinIndexedVector *regionSparse2, CoinIndexedVector *regionSparse3, bool noPermuteRegion3=false)
Updates one column (FTRAN) from region2 Tries to do FT update number returned is negative if no room.
int factorSparseSmall()
Does sparse phase of factorization (for smaller problems) return code is <0 error,...
int messageLevel() const
Level of detail of messages.
double zeroTolerance_
Zero tolerance.
void setCollectStatistics(bool onOff) const
For statistics.
int denseThreshold_
Dense threshold.
double pivotTolerance() const
Pivot tolerance.
int * startRowL() const
Start of each row in L.
void updateColumnTransposeRSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR) when sparse.
void checkConsistency()
Checks that row and column copies look OK.
bool getColumnSpaceIterateR(int iColumn, double value, int iRow)
getColumnSpaceIterateR.
int * startColumnU() const
Start of each column in U.
int maximumRowsExtra_
Maximum number of Rows after iterating.
int sparseThreshold() const
get sparse threshold
int lengthAreaR_
length of area reserved for R
int updateColumn(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2, bool noPermute=false) const
This version has same effect as above with FTUpdate==false so number returned is always >=0.
int biasLU() const
L to U bias 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias.
int factorDense()
Does dense phase of factorization return code is <0 error, 0= finished.
void setDenseThreshold(int value)
Sets dense threshold.
double maximumCoefficient() const
Returns maximum absolute value in factorization.
int numberFtranCounts_
We can roll over factorizations.
int numberColumns_
Number of Columns in factorization.
int deleteColumn(int Row)
Deletes one Column from basis, returns rank.
double conditionNumber() const
Condition number - product of pivots after factorization.
int status_
Status of factorization.
int factorElements_
Number of elements after factorization.
CoinFactorization(const CoinFactorization &other)
Copy constructor.
int * lastRow() const
Returns address of lastRow region.
void gutsOfCopy(const CoinFactorization &other)
CoinIntArrayWithLength nextCount_
Next Row/Column with count.
int deleteRow(int Row)
Deletes one Row from basis, returns rank.
int factorizePart1(int numberRows, int numberColumns, int estimateNumberElements, int *indicesRow[], int *indicesColumn[], CoinFactorizationDouble *elements[], double areaFactor=0.0)
Two part version for maximum flexibility This part creates arrays for user to fill.
void setForrestTomlin(bool value)
int biggerDimension_
Larger of row and column size.
void emptyRows(int numberToEmpty, const int which[])
Takes out all entries for given rows.
int updateColumnUDensish(double *COIN_RESTRICT region, int *COIN_RESTRICT regionIndex) const
Updates part of column (FTRANU)
int denseThreshold() const
Gets dense threshold.
void sparseThreshold(int value)
set sparse threshold
int numberRows_
Number of Rows in factorization.
int numberElementsR() const
Returns number in R area.
int replaceRow(int whichRow, int numberElements, const int indicesColumn[], const double elements[])
Replaces one Row in basis, At present assumes just a singleton on row is in basis returns 0=OK,...
void updateColumnTransposeLDensish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when densish by column.
CoinIntArrayWithLength lastCount_
Previous Row/Column with count.
void updateOneColumnTranspose(CoinIndexedVector *regionWork, int &statistics) const
Part of twocolumnsTranspose.
int * densePermute_
Dense permutation.
CoinIntArrayWithLength permute_
Permutation vector for pivot row order.
bool getColumnSpace(int iColumn, int extraNeeded)
Gets space for one Column with given length, may have to do compression (returns True if successful),...
void updateColumnTransposeL(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL)
int * indexRowL() const
Row indices of L.
void updateColumnTransposeUByColumn(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) by column assumes index is sorted i.e.
int saveFactorization(const char *file) const
Debug - save on file - 0 if no error.
void sort() const
Debug - sort so can compare.
int * numberInColumn() const
Number of entries in each column.
CoinFactorizationDoubleArrayWithLength pivotRegion_
Inverses of pivot values.
CoinIntArrayWithLength indexRowU_
Row indices of U.
int numberGoodColumns() const
Number of good columns in factorization.
int numberR_
Number in R.
CoinIntArrayWithLength nextRow_
Next Row in memory order.
int maximumU_
Maximum space used in U.
int * permute() const
Returns address of permute region.
int factorSparseLarge()
Does sparse phase of factorization (for larger problems) return code is <0 error, 0= finished.
int sparseThreshold2_
And one for "sparsish".
double slackValue_
Whether slack value is +1 or -1.
bool pivotRowSingleton(int pivotRow, int pivotColumn)
Does one pivot on Row Singleton in factorization.
int numberTrials_
0 - no increasing rows - no permutations, 1 - no increasing rows but permutations 2 - increasing rows
void deleteLink(int index)
Deletes a link in chain of equal counts.
double pivotTolerance_
Pivot tolerance.
void setPivots(int value)
Sets number of pivots since factorization.
CoinFactorizationDouble * elementU() const
Elements of U.
int lengthAreaL() const
Returns length of L area.
void updateColumnRFT(CoinIndexedVector *region, int *indexIn)
Updates part of column (FTRANR) with FT update.
void areaFactor(double value)
CoinIntArrayWithLength convertRowToColumnU_
Converts rows to columns in U.
CoinFactorization & operator=(const CoinFactorization &other)
= copy
double slackValue() const
Whether slack value is +1 or -1.
void setStatus(int value)
Sets status.
CoinIntArrayWithLength numberInRow_
Number in each Row.
int * indexColumnL() const
Index of column in row for L.
int lengthL_
Length of L.
CoinIntArrayWithLength indexColumnU_
Base address for U (may change)
void almostDestructor()
Delete all stuff (leaves as after CoinFactorization())
void updateColumnPFI(CoinIndexedVector *regionSparse) const
Updates part of column PFI (FTRAN) (after rest)
CoinIntArrayWithLength startColumnU_
Start of each column in U.
int numberColumns() const
Total number of columns in factorization.
void updateColumnLSparsish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when sparsish.
void gutsOfDestructor(int type=1)
The real work of constructors etc 0 just scalars, 1 bit normal.
int numberDense() const
Returns number of dense rows.
CoinIntArrayWithLength firstCount_
First Row/Column with count of k, can tell which by offset - Rows then Columns.
int biasLU_
L to U bias 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias.
int factorize(int numberRows, int numberColumns, int numberElements, int maximumL, int maximumU, const int indicesRow[], const int indicesColumn[], const double elements[], int permutation[], double areaFactor=0.0)
When given as triplets.
void updateColumnUSparse(CoinIndexedVector *regionSparse, int *indexIn) const
Updates part of column (FTRANU) when sparse.
void replaceColumnU(CoinIndexedVector *regionSparse, int *deleted, int internalPivotRow)
Combines BtranU and delete elements If deleted is NULL then delete elements otherwise store where ele...
bool pivotOneOtherRow(int pivotRow, int pivotColumn)
Pivots when just one other row so faster?
int numberElementsU() const
Returns number in U area.
int * pivotColumnBack() const
Returns address of pivotColumnBack region (also used for permuting) Now uses firstCount to save memor...
int numberRowsExtra() const
Number of Rows after iterating.
CoinIntArrayWithLength permuteBack_
DePermutation vector for pivot row order.
int addColumn(int numberElements, int indicesRow[], double elements[])
Adds one Column to basis, can increase size of basis.
double * denseAreaAddress_
Dense area - actually used (for alignment etc)
int numberSlacks_
Number of slacks at beginning of U.
CoinIntArrayWithLength nextColumn_
Next Column in memory order.
void slackValue(double value)
int lengthR_
Length of R stuff.
int numberForrestTomlin() const
Length of FT vector.
int numberPivots_
Number pivots since last factorization.
void goSparse()
makes a row copy of L for speed and to allow very sparse problems
CoinIntArrayWithLength numberInColumnPlus_
Number in each Column including pivoted.
int factorSparse()
Does sparse phase of factorization return code is <0 error, 0= finished.
void resetStatistics()
Reset all sparsity etc statistics.
bool collectStatistics() const
For statistics.
double relaxCheck_
Relax check on accuracy in replaceColumn.
int add(int numberElements, int indicesRow[], int indicesColumn[], double elements[])
Adds given elements to Basis and updates factorization, can increase size of basis.
int * startColumnL() const
Start of each column in L.
int numberL() const
Number in L.
void setPersistenceFlag(int value)
bool doForrestTomlin_
true if Forrest Tomlin update, false if PFI
void setNumberRows(int value)
Set number of Rows after factorization.
int lengthU_
Base of U is always 0.
void preProcess(int state, int possibleDuplicates=-1)
PreProcesses raw triplet data.
int * numberInRow() const
Number of entries in each row.
double zeroTolerance() const
Zero tolerance.
bool pivot(int pivotRow, int pivotColumn, int pivotRowPosition, int pivotColumnPosition, CoinFactorizationDouble work[], unsigned int workArea2[], int increment2, T markRow[], int largeInteger)
int pivots() const
Returns number of pivots since factorization.
void updateColumnTransposeRDensish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR) when dense.
CoinIntArrayWithLength lastColumn_
Previous Column in memory order.
void addLink(int index, int count)
Adds a link in chain of equal counts.
int factorize(const CoinPackedMatrix &matrix, int rowIsBasic[], int columnIsBasic[], double areaFactor=0.0)
When part of LP - given by basic variables.
int addRow(int numberElements, int indicesColumn[], double elements[])
Adds one Row to basis, can increase size of basis.
CoinIntArrayWithLength startColumnR_
Start of columns for R.
CoinFactorizationDoubleArrayWithLength elementL_
Elements of L.
void updateColumnTransposePFI(CoinIndexedVector *region) const
Updates part of column transpose PFI (BTRAN) (before rest)
CoinIntArrayWithLength lastRow_
Previous Row in memory order.
bool pivotColumnSingleton(int pivotRow, int pivotColumn)
Does one pivot on Column Singleton in factorization.
CoinFactorizationDoubleArrayWithLength workArea_
First work area.
void updateTwoColumnsTranspose(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2, CoinIndexedVector *regionSparse3, int type) const
Updates two columns (BTRAN) from regionSparse2 and 3 regionSparse starts as zero and is zero at end N...
int lengthAreaU() const
Returns length of U area.
int persistenceFlag() const
Array persistence flag If 0 then as now (delete/new) 1 then only do arrays if bigger needed 2 as 1 bu...
int persistenceFlag_
Array persistence flag If 0 then as now (delete/new) 1 then only do arrays if bigger needed 2 as 1 bu...
CoinIntArrayWithLength indexRowL_
Row indices of L.
void clearArrays()
Get rid of all memory.
void setBiasLU(int value)
int numberU_
Number in U.
int totalElements_
Number of elements in U (to go) or while iterating total overall.
void gutsOfInitialize(int type)
1 bit - tolerances etc, 2 more, 4 dummy arrays
void updateColumnR(CoinIndexedVector *region) const
Updates part of column (FTRANR) without FT update.
double areaFactor() const
Whether larger areas needed.
int factorizePart2(int permutation[], int exactNumberElements)
This is part two of factorization Arrays belong to factorization and were returned by part 1 If statu...
int * indexRowR_
Row indices for R.
int sparseThreshold_
Below this use sparse technology - if 0 then no L row copy.
CoinFactorizationDouble * pivotRegion() const
Returns address of pivot region.
int replaceColumnPFI(CoinIndexedVector *regionSparse, int pivotRow, double alpha)
Replaces one Column to basis for PFI returns 0=OK, 1=Probably OK, 2=singular, 3=no room.
int numberCompressions() const
Number of compressions done.
bool getRowSpace(int iRow, int extraNeeded)
Gets space for one Row with given length, may have to do compression (returns True if successful),...
int maximumPivots() const
Maximum number of pivots between factorizations.
int factor()
Does most of factorization.
int updateColumnFT(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2)
Updates one column (FTRAN) from regionSparse2 Tries to do FT update number returned is negative if no...
CoinIntArrayWithLength startRowU_
Start of each Row as pointer.
int * array() const
Get Array.
Sparse Matrix Base Class.