Cbc  2.10.5
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CoinFactorization.hpp
Go to the documentation of this file.
1 /* $Id: CoinFactorization.hpp 2083 2019-01-06 19:38:09Z unxusr $ */
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 
25 class CoinPackedMatrix;
51  friend void CoinFactorizationUnitTest(const std::string &mpsDir);
52 
53 public:
60 
64  void almostDestructor();
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);
101  int factorize(int numberRows,
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);
114  int factorizePart1(int numberRows,
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 
135  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
324  inline double slackValue() const
326  {
327  return slackValue_;
328  }
329  void slackValue(double value);
330 #endif
331  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  }
421  inline int maximumColumnsExtra()
422  {
423  return maximumColumnsExtra_;
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 
458  int replaceColumn(CoinIndexedVector *regionSparse,
459  int pivotRow,
460  double pivotCheck,
461  bool checkBeforeModifying = false,
462  double acceptablePivot = 1.0e-8);
467  void replaceColumnU(CoinIndexedVector *regionSparse,
468  int *deleted,
469  int internalPivotRow);
470 #ifdef ABC_USE_COIN_FACTORIZATION
471 
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
580  void replaceColumn1(CoinIndexedVector *regionSparse, int pivotRow);
583  int replaceColumn2(CoinIndexedVector *regionSparse,
584  int pivotRow,
585  double pivotCheck);
586 #endif
587 
588 
598  int updateColumnFT(CoinIndexedVector *regionSparse,
599  CoinIndexedVector *regionSparse2);
602  int updateColumn(CoinIndexedVector *regionSparse,
603  CoinIndexedVector *regionSparse2,
604  bool noPermute = false) const;
610  int updateTwoColumnsFT(CoinIndexedVector *regionSparse1,
611  CoinIndexedVector *regionSparse2,
612  CoinIndexedVector *regionSparse3,
613  bool noPermuteRegion3 = false);
618  int updateColumnTranspose(CoinIndexedVector *regionSparse,
619  CoinIndexedVector *regionSparse2) const;
621  void updateOneColumnTranspose(CoinIndexedVector *regionWork, int &statistics) const;
626  void updateTwoColumnsTranspose(CoinIndexedVector *regionSparse,
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 
645  inline void clearArrays()
647  {
649  }
651 
654 
657  int add(int numberElements,
658  int indicesRow[],
659  int indicesColumn[], double elements[]);
660 
663  int addColumn(int numberElements,
664  int indicesRow[], double elements[]);
665 
668  int addRow(int numberElements,
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 
685  void checkSparse();
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
704  void gutsOfDestructor(int type = 1);
707  void gutsOfInitialize(int type);
708  void gutsOfCopy(const CoinFactorization &other);
709 
711  void resetStatistics();
712 
714 
716  void getAreas(int numberRows,
718  int numberColumns,
719  int maximumL,
720  int maximumU);
721 
724  void preProcess(int state,
725  int possibleDuplicates = -1);
727  int factor();
728 
729 protected:
732  int factorSparse();
735  int factorSparseSmall();
738  int factorSparseLarge();
741  int factorDense();
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);
786  void checkConsistency();
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 
851  void updateColumnUSparse(CoinIndexedVector *regionSparse,
852  int *indexIn) const;
854  void updateColumnUSparsish(CoinIndexedVector *regionSparse,
855  int *indexIn) const;
857  int updateColumnUDensish(double *COIN_RESTRICT region,
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 
874  void updateColumnTransposePFI(CoinIndexedVector *region) const;
878  int smallestIndex) const;
882  int smallestIndex) const;
886  int smallestIndex) const;
893  int smallestIndex) const;
894 
896  void updateColumnTransposeR(CoinIndexedVector *region) const;
901 
903  void updateColumnTransposeL(CoinIndexedVector *region) const;
912 
913 public:
918  int replaceColumnPFI(CoinIndexedVector *regionSparse,
919  int pivotRow, double alpha);
920 
921 protected:
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();
953  int *numberInRow = numberInRow_.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 
1010  int *startColumnL = startColumnL_.array();
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)
1103  numberInColumn[pivotColumn] = 0;
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);
1408  deleteLink(pivotColumn + numberRows_);
1409  totalElements_ += added;
1410  return true;
1411  }
1412 
1413  /********************************* END LARGE TEMPLATE ********/
1415 protected:
1419  double pivotTolerance_;
1423 #ifndef COIN_FAST_CODE
1424  double slackValue_;
1426 #else
1427 #ifndef slackValue_
1428 #define slackValue_ -1.0
1429 #endif
1430 #endif
1431  double areaFactor_;
1434  double relaxCheck_;
1469  int status_;
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 
1568  int baseL_;
1569 
1572 
1575 
1578 
1581 
1584 
1587 
1590 
1593 
1596 
1599 
1602 
1605 
1607  double *denseArea_;
1608 
1611 
1614 
1617 
1620 
1623 
1626 
1629 
1630 public:
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 
1642  mutable int numberFtranCounts_;
1643  mutable int numberBtranCounts_;
1644 
1652 
1653 protected:
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
1679  int sparseOffset_;
1681 #endif
1682 
1685  int biasLU_;
1692 #ifdef ABC_USE_COIN_FACTORIZATION
1693  int parallelMode_;
1695 #endif
1696 
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 ;
1711 typedef int ipfint;
1712 typedef 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 
1769  int *startColumnL = startColumnL_.array();
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  }
2170 BAD_PIVOT:
2171 
2172  ;
2173 }
2174 #undef FAC_UNSET
2175 #endif
2176 
2177 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2178 */
CoinIntArrayWithLength startRowL_
Start of each row in L.
int numberR_
Number in R.
#define COINFACTORIZATION_BITS_PER_INT
int numberCompressions_
Number of compressions done.
CoinFactorizationDouble * elementU() const
Elements of U.
bool getColumnSpaceIterateR(int iColumn, double value, int iRow)
getColumnSpaceIterateR.
CoinIntArrayWithLength nextColumn_
Next Column in memory order.
int lengthAreaL() const
Returns length of L area.
bool pivotOneOtherRow(int pivotRow, int pivotColumn)
Pivots when just one other row so faster?
CoinFactorizationDoubleArrayWithLength pivotRegion_
Inverses of pivot values.
void setPivots(int value)
Sets number of pivots since factorization.
void separateLinks(int count, bool rowsFirst)
Separate out links with same row/column count.
int totalElements_
Number of elements in U (to go) or while iterating total overall.
bool doForrestTomlin_
true if Forrest Tomlin update, false if PFI
CoinFactorizationDouble * pivotRegion() const
Returns address of pivot region.
void updateColumnTransposeLDensish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when densish by column.
int numberFtranCounts_
We can roll over factorizations.
int numberPivots_
Number pivots since last factorization.
void updateColumnTransposeRSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR) when sparse.
double areaFactor() const
Whether larger areas needed.
int numberDense() const
Returns number of dense rows.
int maximumRowsExtra_
Maximum number of Rows after iterating.
int lengthL_
Length of L.
int status() const
Returns status.
int baseL() const
Base of L.
int denseThreshold() const
Gets dense threshold.
int maximumColumnsExtra()
Maximum number of Columns after iterating.
int sparseThreshold_
Below this use sparse technology - if 0 then no L row copy.
int * indexColumnL() const
Index of column in row for L.
int getColumnSpaceIterate(int iColumn, double value, int iRow)
getColumnSpaceIterate.
void deleteLink(int index)
Deletes a link in chain of equal counts.
int numberElementsL() const
Returns number in L area.
void checkSparse()
See if worth going sparse.
int * pivotColumn() const
Returns address of pivotColumn region (also used for permuting)
CoinIntArrayWithLength lastColumn_
Previous Column in memory order.
CoinFactorizationDoubleArrayWithLength elementL_
Elements of L.
void setNumberRows(int value)
Set number of Rows after factorization.
bool collectStatistics() const
For statistics.
int lengthAreaL_
Length of area reserved for L.
void setForrestTomlin(bool value)
void updateColumnTransposeUSparsish(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) when sparsish, assumes index is sorted i.e.
int status_
Status of factorization.
#define COINFACTORIZATION_SHIFT_PER_INT
void updateColumnTransposeR(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR)
CoinIntArrayWithLength pivotRowL_
Pivots for L.
int * startColumnL() const
Start of each column in L.
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...
bool pivotRowSingleton(int pivotRow, int pivotColumn)
Does one pivot on Row Singleton in factorization.
void relaxAccuracyCheck(double value)
Allows change of pivot accuracy check 1.0 == none &gt;1.0 relaxed.
void sort() const
Debug - sort so can compare.
int numberRows() const
Number of Rows after factorization.
~CoinFactorization()
Destructor.
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 &gt;=0...
double maximumCoefficient() const
Returns maximum absolute value in factorization.
int * permute() const
Returns address of permute region.
int numberGoodU_
Number factorized in U (not row singletons)
int * startColumnU() const
Start of each column in U.
CoinFactorization()
Default constructor.
int add(int numberElements, int indicesRow[], int indicesColumn[], double elements[])
Adds given elements to Basis and updates factorization, can increase size of basis.
int deleteRow(int Row)
Deletes one Row from basis, returns rank.
double * denseArea_
Dense area.
CoinIntArrayWithLength markRow_
Marks rows to be updated.
int * startRowL() const
Start of each row in L.
int numberRowsExtra() const
Number of Rows after iterating.
void gutsOfInitialize(int type)
1 bit - tolerances etc, 2 more, 4 dummy arrays
int numberElements() const
Total number of elements in factorization.
void updateColumnTransposeLByRow(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when densish by row.
int updateColumnTranspose(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2) const
Updates one column (BTRAN) from regionSparse2 regionSparse starts as zero and is zero at end Note - i...
void replaceColumnU(CoinIndexedVector *regionSparse, int *deleted, int internalPivotRow)
Combines BtranU and delete elements If deleted is NULL then delete elements otherwise store where ele...
CoinIntArrayWithLength startRowU_
Start of each Row as pointer.
int updateColumnFT(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2)
Updates one column (FTRAN) from regionSparse2 Tries to do FT update number returned is negative if no...
void updateColumnL(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL)
int numberColumns() const
Total number of columns in factorization.
int numberL_
Number in L.
double relaxCheck_
Relax check on accuracy in replaceColumn.
int numberL() const
Number in L.
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.
void cleanup()
Cleans up at end of factorization.
int addColumn(int numberElements, int indicesRow[], double elements[])
Adds one Column to basis, can increase size of basis.
bool getRowSpaceIterate(int iRow, int extraNeeded)
Gets space for one Row with given length while iterating, may have to do compression (returns True i...
int * numberInRow() const
Number of entries in each row.
int restoreFactorization(const char *file, bool factor=false)
Debug - restore from file - 0 if no error on file.
CoinFactorization & operator=(const CoinFactorization &other)
= copy
int lengthAreaU() const
Returns length of U area.
int * pivotColumnBack() const
Returns address of pivotColumnBack region (also used for permuting) Now uses firstCount to save memor...
void updateColumnTransposeRDensish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR) when dense.
int checkPivot(double saveFromU, double oldPivot) const
Returns accuracy status of replaceColumn returns 0=OK, 1=Probably OK, 2=singular. ...
double zeroTolerance() const
Zero tolerance.
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 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 updateColumnTransposeLSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparse (by Row)
#define COIN_RESTRICT
double * denseAreaAddress_
Dense area - actually used (for alignment etc)
void updateOneColumnTranspose(CoinIndexedVector *regionWork, int &statistics) const
Part of twocolumnsTranspose.
void updateColumnPFI(CoinIndexedVector *regionSparse) const
Updates part of column PFI (FTRAN) (after rest)
double pivotTolerance_
Pivot tolerance.
#define collectStatistics_
For statistics.
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 saveFactorization(const char *file) const
Debug - save on file - 0 if no error.
void gutsOfCopy(const CoinFactorization &other)
int * numberInColumn() const
Number of entries in each column.
int factorSparse()
Does sparse phase of factorization return code is &lt;0 error, 0= finished.
int numberRowsExtra_
Number of Rows after iterating.
CoinIntArrayWithLength numberInColumn_
Number in each Column.
void updateColumnLSparse(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when sparse.
int sparseThreshold() const
get sparse threshold
int persistenceFlag_
Array persistence flag If 0 then as now (delete/new) 1 then only do arrays if bigger needed 2 as 1 bu...
int factor()
Does most of factorization.
void setPersistenceFlag(int value)
int numberGoodColumns() const
Number of good columns in factorization.
CoinFactorizationDouble * array() const
Get Array.
int addRow(int numberElements, int indicesColumn[], double elements[])
Adds one Row to basis, can increase size of basis.
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...
void getAreas(int numberRows, int numberColumns, int maximumL, int maximumU)
Gets space for a factorization, called by constructors.
double CoinFactorizationDouble
Definition: CoinTypes.hpp:57
int numberDense_
Number of dense rows.
void updateColumnLSparsish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when sparsish.
void addLink(int index, int count)
Adds a link in chain of equal counts.
void gutsOfDestructor(int type=1)
The real work of constructors etc 0 just scalars, 1 bit normal.
int deleteColumn(int Row)
Deletes one Column from basis, returns rank.
int * densePermute_
Dense permutation.
void updateColumnTransposePFI(CoinIndexedVector *region) const
Updates part of column transpose PFI (BTRAN) (before rest)
int messageLevel() const
Level of detail of messages.
double getAccuracyCheck() const
void setBiasLU(int value)
int pivots() const
Returns number of pivots since factorization.
Indexed Vector.
CoinIntArrayWithLength nextCount_
Next Row/Column with count.
bool pivotColumnSingleton(int pivotRow, int pivotColumn)
Does one pivot on Column Singleton in factorization.
void updateColumnTransposeLSparsish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparsish by row.
CoinIntArrayWithLength lastRow_
Previous Row in memory order.
int denseThreshold_
Dense threshold.
void almostDestructor()
Delete all stuff (leaves as after CoinFactorization())
friend void CoinFactorizationUnitTest(const std::string &mpsDir)
int lengthAreaU_
Length of area reserved for U.
#define COINFACTORIZATION_MASK_PER_INT
CoinIntArrayWithLength firstCount_
First Row/Column with count of k, can tell which by offset - Rows then Columns.
int * lastRow() const
Returns address of lastRow region.
CoinIntArrayWithLength sparse_
Sparse regions.
CoinFactorizationDouble * elementR_
Elements of R.
void updateColumnTransposeL(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL)
void updateColumnRFT(CoinIndexedVector *region, int *indexIn)
Updates part of column (FTRANR) with FT update.
CoinIntArrayWithLength startColumnL_
Start of each column in L.
int * indexRowL() const
Row indices of L.
void setStatus(int value)
Sets status.
double ftranAverageAfterL_
While these are average ratios collected over last period.
CoinIntArrayWithLength pivotColumn_
Pivot order for each Column.
void show_self() const
Debug show object (shows one representation)
double zeroTolerance_
Zero tolerance.
int numberU_
Number in U.
int biggerDimension_
Larger of row and column size.
int biasLU_
L to U bias 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias.
void resetStatistics()
Reset all sparsity etc statistics.
void areaFactor(double value)
int maximumPivots() const
Maximum number of pivots between factorizations.
double adjustedAreaFactor() const
Returns areaFactor but adjusted for dense.
int numberElementsR() const
Returns number in R area.
int factorSparseLarge()
Does sparse phase of factorization (for larger problems) return code is &lt;0 error, 0= finished...
double ftranCountInput_
Below are all to collect.
int numberCompressions() const
Number of compressions done.
int numberTrials_
0 - no increasing rows - no permutations, 1 - no increasing rows but permutations 2 - increasing rows...
bool pivot(int pivotRow, int pivotColumn, int pivotRowPosition, int pivotColumnPosition, CoinFactorizationDouble work[], unsigned int workArea2[], int increment2, T markRow[], int largeInteger)
void updateColumnUSparsish(CoinIndexedVector *regionSparse, int *indexIn) const
Updates part of column (FTRANU) when sparsish.
CoinIntArrayWithLength nextRow_
Next Row in memory order.
CoinFactorizationDoubleArrayWithLength elementU_
Elements of U.
double slackValue() const
Whether slack value is +1 or -1.
Sparse Matrix Base Class.
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 numberColumns_
Number of Columns in factorization.
CoinIntArrayWithLength indexColumnU_
Base address for U (may change)
int lengthR_
Length of R stuff.
void updateColumnUSparse(CoinIndexedVector *regionSparse, int *indexIn) const
Updates part of column (FTRANU) when sparse.
CoinFactorizationDoubleArrayWithLength workArea_
First work area.
CoinIntArrayWithLength saveColumn_
Columns left to do in a single pivot.
CoinIntArrayWithLength convertRowToColumnU_
Converts rows to columns in U.
int * indexRowR_
Row indices for R.
void clearArrays()
Get rid of all memory.
void updateColumnTransposeUDensish(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) when densish, assumes index is sorted i.e.
int factorDense()
Does dense phase of factorization return code is &lt;0 error, 0= finished.
CoinIntArrayWithLength numberInColumnPlus_
Number in each Column including pivoted.
int updateColumnUDensish(double *COIN_RESTRICT region, int *COIN_RESTRICT regionIndex) const
Updates part of column (FTRANU)
CoinIntArrayWithLength pivotColumnBack_
Inverse Pivot order for each Column.
CoinIntArrayWithLength permuteBack_
DePermutation vector for pivot row order.
double slackValue_
Whether slack value is +1 or -1.
void goSparse()
makes a row copy of L for speed and to allow very sparse problems
CoinFactorizationDouble * version.
CoinIntArrayWithLength indexColumnL_
Index of column in row for L.
CoinIntArrayWithLength indexRowL_
Row indices of L.
CoinIntArrayWithLength permute_
Permutation vector for pivot row order.
void updateColumnTransposeUSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANU) when sparse, assumes index is sorted i.e.
int * indexRowU() const
Row indices of U.
CoinUnsignedIntArrayWithLength workArea2_
Second work area.
This deals with Factorization and Updates.
bool getColumnSpace(int iColumn, int extraNeeded)
Gets space for one Column with given length, may have to do compression (returns True if successful)...
void updateColumnR(CoinIndexedVector *region) const
Updates part of column (FTRANR) without FT update.
double areaFactor_
How much to multiply areas by.
void checkConsistency()
Checks that row and column copies look OK.
int * array() const
Get Array.
void updateColumnU(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANU)
void updateColumnTransposeU(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU), assumes index is sorted i.e.
int numberSlacks_
Number of slacks at beginning of U.
int numberColumnsExtra_
Number of Columns after iterating.
int factorElements_
Number of elements after factorization.
CoinIntArrayWithLength numberInRow_
Number in each Row.
void setNumberElementsU(int value)
Setss number in U area.
int lengthAreaR_
length of area reserved for R
int numberForrestTomlin() const
Length of FT vector.
void setDenseThreshold(int value)
Sets dense threshold.
int messageLevel_
Detail in messages.
CoinIntArrayWithLength lastCount_
Previous Row/Column with count.
void updateColumnTransposeUByColumn(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) by column assumes index is sorted i.e.
CoinIntArrayWithLength startColumnR_
Start of columns for R.
CoinFactorizationDouble * elementByRowL() const
Elements in L (row copy)
int maximumPivots_
Maximum number of pivots before factorization.
bool spaceForForrestTomlin() const
True if FT update and space.
bool reorderU()
Reorders U so contiguous and in order (if there is space) Returns true if it could.
void preProcess(int state, int possibleDuplicates=-1)
PreProcesses raw triplet data.
int maximumColumnsExtra_
Maximum number of Columns after iterating.
int maximumRowsExtra() const
Maximum of Rows after iterating.
int numberElementsU() const
Returns number in 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 biasLU() const
L to U bias 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias.
int sparseThreshold2_
And one for &quot;sparsish&quot;.
int numberGoodL_
Number factorized in L.
double conditionNumber() const
Condition number - product of pivots after factorization.
double pivotTolerance() const
Pivot tolerance.
CoinFactorizationDoubleArrayWithLength elementByRowL_
Elements in L (row copy)
int numberRows_
Number of Rows in factorization.
void setCollectStatistics(bool onOff) const
For statistics.
CoinIntArrayWithLength startColumnU_
Start of each column in U.
void updateColumnLDensish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when densish.
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...
int factorSparseSmall()
Does sparse phase of factorization (for smaller problems) return code is &lt;0 error, 0= finished.
CoinIntArrayWithLength indexRowU_
Row indices of U.
bool forrestTomlin() const
true if Forrest Tomlin update, false if PFI
void emptyRows(int numberToEmpty, const int which[])
Takes out all entries for given rows.
int factorize(const CoinPackedMatrix &matrix, int rowIsBasic[], int columnIsBasic[], double areaFactor=0.0)
When part of LP - given by basic variables.
bool getRowSpace(int iRow, int extraNeeded)
Gets space for one Row with given length, may have to do compression (returns True if successful)...
int * permuteBack() const
Returns address of permuteBack region.
int maximumU_
Maximum space used in U.
int lengthU_
Base of U is always 0.