CoinFactorization.hpp
Go to the documentation of this file.
1 /* $Id: CoinFactorization.hpp 1277 2010-04-21 20:19:32Z forrest $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 
5 /*
6  Authors
7 
8  John Forrest
9 
10  */
11 #ifndef CoinFactorization_H
12 #define CoinFactorization_H
13 //#define COIN_ONE_ETA_COPY 100
14 
15 #include <iostream>
16 #include <string>
17 #include <cassert>
18 #include <cstdio>
19 #include "CoinFinite.hpp"
20 #include "CoinIndexedVector.hpp"
21 class CoinPackedMatrix;
48  friend void CoinFactorizationUnitTest( const std::string & mpsDir );
49 
50 public:
51 
57  CoinFactorization ( const CoinFactorization &other);
58 
62  void almostDestructor();
64  void show_self ( ) const;
66  int saveFactorization (const char * file ) const;
70  int restoreFactorization (const char * file , bool factor=false) ;
72  void sort ( ) const;
76 
86  int factorize ( const CoinPackedMatrix & matrix,
87  int rowIsBasic[], int columnIsBasic[] ,
88  double areaFactor = 0.0 );
99  int factorize ( int numberRows,
100  int numberColumns,
102  CoinBigIndex maximumL,
103  CoinBigIndex maximumU,
104  const int indicesRow[],
105  const int indicesColumn[], const double elements[] ,
106  int permutation[],
107  double areaFactor = 0.0);
112  int factorizePart1 ( int numberRows,
113  int numberColumns,
114  CoinBigIndex estimateNumberElements,
115  int * indicesRow[],
116  int * indicesColumn[],
117  CoinFactorizationDouble * elements[],
118  double areaFactor = 0.0);
125  int factorizePart2 (int permutation[],int exactNumberElements);
127  double conditionNumber() const;
128 
130 
133  inline int status ( ) const {
135  return status_;
136  }
138  inline void setStatus ( int value)
139  { status_=value; }
141  inline int pivots ( ) const {
142  return numberPivots_;
143  }
145  inline void setPivots ( int value )
146  { numberPivots_=value; }
148  inline int *permute ( ) const {
149  return permute_.array();
150  }
152  inline int *pivotColumn ( ) const {
153  return pivotColumn_.array();
154  }
157  return pivotRegion_.array();
158  }
160  inline int *permuteBack ( ) const {
161  return permuteBack_.array();
162  }
165  inline int *pivotColumnBack ( ) const {
166  //return firstCount_.array();
167  return pivotColumnBack_.array();
168  }
170  inline CoinBigIndex * startRowL() const
171  { return startRowL_.array();}
172 
174  inline CoinBigIndex * startColumnL() const
175  { return startColumnL_.array();}
176 
178  inline int * indexColumnL() const
179  { return indexColumnL_.array();}
180 
182  inline int * indexRowL() const
183  { return indexRowL_.array();}
184 
187  { return elementByRowL_.array();}
188 
190  inline int numberRowsExtra ( ) const {
191  return numberRowsExtra_;
192  }
194  inline void setNumberRows(int value)
195  { numberRows_ = value; }
197  inline int numberRows ( ) const {
198  return numberRows_;
199  }
201  inline CoinBigIndex numberL() const
202  { return numberL_;}
203 
205  inline CoinBigIndex baseL() const
206  { return baseL_;}
208  inline int maximumRowsExtra ( ) const {
209  return maximumRowsExtra_;
210  }
212  inline int numberColumns ( ) const {
213  return numberColumns_;
214  }
216  inline int numberElements ( ) const {
217  return totalElements_;
218  }
220  inline int numberForrestTomlin ( ) const {
222  }
224  inline int numberGoodColumns ( ) const {
225  return numberGoodU_;
226  }
228  inline double areaFactor ( ) const {
229  return areaFactor_;
230  }
231  inline void areaFactor ( double value ) {
232  areaFactor_=value;
233  }
235  double adjustedAreaFactor() const;
237  inline void relaxAccuracyCheck(double value)
238  { relaxCheck_ = value;}
239  inline double getAccuracyCheck() const
240  { return relaxCheck_;}
242  inline int messageLevel ( ) const {
243  return messageLevel_ ;
244  }
245  void messageLevel ( int value );
247  inline int maximumPivots ( ) const {
248  return maximumPivots_ ;
249  }
250  void maximumPivots ( int value );
251 
253  inline int denseThreshold() const
254  { return denseThreshold_;}
256  inline void setDenseThreshold(int value)
257  { denseThreshold_ = value;}
259  inline double pivotTolerance ( ) const {
260  return pivotTolerance_ ;
261  }
262  void pivotTolerance ( double value );
264  inline double zeroTolerance ( ) const {
265  return zeroTolerance_ ;
266  }
267  void zeroTolerance ( double value );
268 #ifndef COIN_FAST_CODE
269  inline double slackValue ( ) const {
271  return slackValue_ ;
272  }
273  void slackValue ( double value );
274 #endif
275  double maximumCoefficient() const;
278  inline bool forrestTomlin() const
279  { return doForrestTomlin_;}
280  inline void setForrestTomlin(bool value)
281  { doForrestTomlin_=value;}
283  inline bool spaceForForrestTomlin() const
284  {
286  CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
287  return (space>=0)&&doForrestTomlin_;
288  }
290 
293 
295  inline int numberDense() const
296  { return numberDense_;}
297 
299  inline CoinBigIndex numberElementsU ( ) const {
300  return lengthU_;
301  }
303  inline void setNumberElementsU(CoinBigIndex value)
304  { lengthU_ = value; }
306  inline CoinBigIndex lengthAreaU ( ) const {
307  return lengthAreaU_;
308  }
310  inline CoinBigIndex numberElementsL ( ) const {
311  return lengthL_;
312  }
314  inline CoinBigIndex lengthAreaL ( ) const {
315  return lengthAreaL_;
316  }
318  inline CoinBigIndex numberElementsR ( ) const {
319  return lengthR_;
320  }
323  { return numberCompressions_;}
325  inline int * numberInRow() const
326  { return numberInRow_.array();}
328  inline int * numberInColumn() const
329  { return numberInColumn_.array();}
332  { return elementU_.array();}
334  inline int * indexRowU() const
335  { return indexRowU_.array();}
337  inline CoinBigIndex * startColumnU() const
338  { return startColumnU_.array();}
340  inline int maximumColumnsExtra()
341  { return maximumColumnsExtra_;}
345  inline int biasLU() const
346  { return biasLU_;}
347  inline void setBiasLU(int value)
348  { biasLU_=value;}
354  inline int persistenceFlag() const
355  { return persistenceFlag_;}
356  void setPersistenceFlag(int value);
358 
361 
369  int replaceColumn ( CoinIndexedVector * regionSparse,
370  int pivotRow,
371  double pivotCheck ,
372  bool checkBeforeModifying=false,
373  double acceptablePivot=1.0e-8);
378  void replaceColumnU ( CoinIndexedVector * regionSparse,
379  CoinBigIndex * deleted,
380  int internalPivotRow);
382 
392  int updateColumnFT ( CoinIndexedVector * regionSparse,
393  CoinIndexedVector * regionSparse2);
396  int updateColumn ( CoinIndexedVector * regionSparse,
397  CoinIndexedVector * regionSparse2,
398  bool noPermute=false) const;
404  int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
405  CoinIndexedVector * regionSparse2,
406  CoinIndexedVector * regionSparse3,
407  bool noPermuteRegion3=false) ;
412  int updateColumnTranspose ( CoinIndexedVector * regionSparse,
413  CoinIndexedVector * regionSparse2) const;
415  void goSparse();
417  inline int sparseThreshold ( ) const
418  { return sparseThreshold_;}
420  void sparseThreshold ( int value );
422 
427  inline void clearArrays()
429  { gutsOfDestructor();}
431 
434 
437  int add ( CoinBigIndex numberElements,
438  int indicesRow[],
439  int indicesColumn[], double elements[] );
440 
443  int addColumn ( CoinBigIndex numberElements,
444  int indicesRow[], double elements[] );
445 
448  int addRow ( CoinBigIndex numberElements,
449  int indicesColumn[], double elements[] );
450 
452  int deleteColumn ( int Row );
454  int deleteRow ( int Row );
455 
459  int replaceRow ( int whichRow, int numberElements,
460  const int indicesColumn[], const double elements[] );
462  void emptyRows(int numberToEmpty, const int which[]);
464 
465  void checkSparse();
468  inline bool collectStatistics() const
469  { return collectStatistics_;}
471  inline void setCollectStatistics(bool onOff) const
472  { collectStatistics_ = onOff;}
474  void gutsOfDestructor(int type=1);
476  void gutsOfInitialize(int type);
477  void gutsOfCopy(const CoinFactorization &other);
478 
480  void resetStatistics();
481 
482 
484 
486  void getAreas ( int numberRows,
488  int numberColumns,
489  CoinBigIndex maximumL,
490  CoinBigIndex maximumU );
491 
494  void preProcess ( int state,
495  int possibleDuplicates = -1 );
497  int factor ( );
498 protected:
501  int factorSparse ( );
504  int factorSparseSmall ( );
507  int factorSparseLarge ( );
510  int factorDense ( );
511 
513  bool pivotOneOtherRow ( int pivotRow,
514  int pivotColumn );
516  bool pivotRowSingleton ( int pivotRow,
517  int pivotColumn );
519  bool pivotColumnSingleton ( int pivotRow,
520  int pivotColumn );
521 
526  bool getColumnSpace ( int iColumn,
527  int extraNeeded );
528 
531  bool reorderU();
535  bool getColumnSpaceIterateR ( int iColumn, double value,
536  int iRow);
542  CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
543  int iRow);
547  bool getRowSpace ( int iRow, int extraNeeded );
548 
552  bool getRowSpaceIterate ( int iRow,
553  int extraNeeded );
555  void checkConsistency ( );
557  inline void addLink ( int index, int count ) {
558  int *nextCount = nextCount_.array();
559  int *firstCount = firstCount_.array();
560  int *lastCount = lastCount_.array();
561  int next = firstCount[count];
562  lastCount[index] = -2 - count;
563  if ( next < 0 ) {
564  //first with that count
565  firstCount[count] = index;
566  nextCount[index] = -1;
567  } else {
568  firstCount[count] = index;
569  nextCount[index] = next;
570  lastCount[next] = index;
571  }}
573  inline void deleteLink ( int index ) {
574  int *nextCount = nextCount_.array();
575  int *firstCount = firstCount_.array();
576  int *lastCount = lastCount_.array();
577  int next = nextCount[index];
578  int last = lastCount[index];
579  if ( last >= 0 ) {
580  nextCount[last] = next;
581  } else {
582  int count = -last - 2;
583 
584  firstCount[count] = next;
585  }
586  if ( next >= 0 ) {
587  lastCount[next] = last;
588  }
589  nextCount[index] = -2;
590  lastCount[index] = -2;
591  return;
592  }
594  void separateLinks(int count,bool rowsFirst);
596  void cleanup ( );
597 
599  void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
601  void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
603  void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
605  void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
606 
608  void updateColumnR ( CoinIndexedVector * region ) const;
611  void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
612 
614  void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
615 
617  void updateColumnUSparse ( CoinIndexedVector * regionSparse,
618  int * indexIn) const;
620  void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
621  int * indexIn) const;
623  int updateColumnUDensish ( double * COIN_RESTRICT region,
624  int * COIN_RESTRICT regionIndex) const;
627  int & numberNonZero1,
628  double * COIN_RESTRICT region1,
629  int * COIN_RESTRICT index1,
630  int & numberNonZero2,
631  double * COIN_RESTRICT region2,
632  int * COIN_RESTRICT index2) const;
634  void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
636  void permuteBack ( CoinIndexedVector * regionSparse,
637  CoinIndexedVector * outVector) const;
638 
640  void updateColumnTransposePFI ( CoinIndexedVector * region) const;
644  int smallestIndex) const;
648  int smallestIndex) const;
652  int smallestIndex) const;
655  void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
659  int smallestIndex) const;
660 
662  void updateColumnTransposeR ( CoinIndexedVector * region ) const;
664  void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
666  void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
667 
669  void updateColumnTransposeL ( CoinIndexedVector * region ) const;
671  void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
673  void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
675  void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
677  void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
678 public:
683  int replaceColumnPFI ( CoinIndexedVector * regionSparse,
684  int pivotRow, double alpha);
685 protected:
688  int checkPivot(double saveFromU, double oldPivot) const;
689  /********************************* START LARGE TEMPLATE ********/
690 #ifdef INT_IS_8
691 #define COINFACTORIZATION_BITS_PER_INT 64
692 #define COINFACTORIZATION_SHIFT_PER_INT 6
693 #define COINFACTORIZATION_MASK_PER_INT 0x3f
694 #else
695 #define COINFACTORIZATION_BITS_PER_INT 32
696 #define COINFACTORIZATION_SHIFT_PER_INT 5
697 #define COINFACTORIZATION_MASK_PER_INT 0x1f
698 #endif
699  template <class T> inline bool
700  pivot ( int pivotRow,
701  int pivotColumn,
702  CoinBigIndex pivotRowPosition,
703  CoinBigIndex pivotColumnPosition,
705  unsigned int workArea2[],
706  int increment2,
707  T markRow[] ,
708  int largeInteger)
709 {
710  int *indexColumnU = indexColumnU_.array();
714  int *indexRowU = indexRowU_.array();
715  CoinBigIndex *startRowU = startRowU_.array();
716  int *numberInRow = numberInRow_.array();
718  int *indexRowL = indexRowL_.array();
719  int *saveColumn = saveColumn_.array();
720  int *nextRow = nextRow_.array();
721  int *lastRow = lastRow_.array() ;
722 
723  //store pivot columns (so can easily compress)
724  int numberInPivotRow = numberInRow[pivotRow] - 1;
725  CoinBigIndex startColumn = startColumnU[pivotColumn];
726  int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
727  CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
728  int put = 0;
729  CoinBigIndex startRow = startRowU[pivotRow];
730  CoinBigIndex endRow = startRow + numberInPivotRow + 1;
731 
732  if ( pivotColumnPosition < 0 ) {
733  for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
734  int iColumn = indexColumnU[pivotColumnPosition];
735  if ( iColumn != pivotColumn ) {
736  saveColumn[put++] = iColumn;
737  } else {
738  break;
739  }
740  }
741  } else {
742  for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
743  saveColumn[put++] = indexColumnU[i];
744  }
745  }
746  assert (pivotColumnPosition<endRow);
747  assert (indexColumnU[pivotColumnPosition]==pivotColumn);
748  pivotColumnPosition++;
749  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
750  saveColumn[put++] = indexColumnU[pivotColumnPosition];
751  }
752  //take out this bit of indexColumnU
753  int next = nextRow[pivotRow];
754  int last = lastRow[pivotRow];
755 
756  nextRow[last] = next;
757  lastRow[next] = last;
758  nextRow[pivotRow] = numberGoodU_; //use for permute
759  lastRow[pivotRow] = -2;
760  numberInRow[pivotRow] = 0;
761  //store column in L, compress in U and take column out
763 
764  if ( l + numberInPivotColumn > lengthAreaL_ ) {
765  //need more memory
766  if ((messageLevel_&4)!=0)
767  printf("more memory needed in middle of invert\n");
768  return false;
769  }
770  //l+=currentAreaL_->elementByColumn-elementL;
771  CoinBigIndex lSave = l;
772 
774  startColumnL[numberGoodL_] = l; //for luck and first time
775  numberGoodL_++;
776  startColumnL[numberGoodL_] = l + numberInPivotColumn;
777  lengthL_ += numberInPivotColumn;
778  if ( pivotRowPosition < 0 ) {
779  for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
780  int iRow = indexRowU[pivotRowPosition];
781  if ( iRow != pivotRow ) {
782  indexRowL[l] = iRow;
783  elementL[l] = elementU[pivotRowPosition];
784  markRow[iRow] = static_cast<T>(l - lSave);
785  l++;
786  //take out of row list
787  CoinBigIndex start = startRowU[iRow];
788  CoinBigIndex end = start + numberInRow[iRow];
789  CoinBigIndex where = start;
790 
791  while ( indexColumnU[where] != pivotColumn ) {
792  where++;
793  } /* endwhile */
794 #if DEBUG_COIN
795  if ( where >= end ) {
796  abort ( );
797  }
798 #endif
799  indexColumnU[where] = indexColumnU[end - 1];
800  numberInRow[iRow]--;
801  } else {
802  break;
803  }
804  }
805  } else {
806  CoinBigIndex i;
807 
808  for ( i = startColumn; i < pivotRowPosition; i++ ) {
809  int iRow = indexRowU[i];
810 
811  markRow[iRow] = static_cast<T>(l - lSave);
812  indexRowL[l] = iRow;
813  elementL[l] = elementU[i];
814  l++;
815  //take out of row list
816  CoinBigIndex start = startRowU[iRow];
817  CoinBigIndex end = start + numberInRow[iRow];
818  CoinBigIndex where = start;
819 
820  while ( indexColumnU[where] != pivotColumn ) {
821  where++;
822  } /* endwhile */
823 #if DEBUG_COIN
824  if ( where >= end ) {
825  abort ( );
826  }
827 #endif
828  indexColumnU[where] = indexColumnU[end - 1];
829  numberInRow[iRow]--;
830  assert (numberInRow[iRow]>=0);
831  }
832  }
833  assert (pivotRowPosition<endColumn);
834  assert (indexRowU[pivotRowPosition]==pivotRow);
835  CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
836  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
837 
838  pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
839  pivotRowPosition++;
840  for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
841  int iRow = indexRowU[pivotRowPosition];
842 
843  markRow[iRow] = static_cast<T>(l - lSave);
844  indexRowL[l] = iRow;
845  elementL[l] = elementU[pivotRowPosition];
846  l++;
847  //take out of row list
848  CoinBigIndex start = startRowU[iRow];
849  CoinBigIndex end = start + numberInRow[iRow];
850  CoinBigIndex where = start;
851 
852  while ( indexColumnU[where] != pivotColumn ) {
853  where++;
854  } /* endwhile */
855 #if DEBUG_COIN
856  if ( where >= end ) {
857  abort ( );
858  }
859 #endif
860  indexColumnU[where] = indexColumnU[end - 1];
861  numberInRow[iRow]--;
862  assert (numberInRow[iRow]>=0);
863  }
864  markRow[pivotRow] = static_cast<T>(largeInteger);
865  //compress pivot column (move pivot to front including saved)
866  numberInColumn[pivotColumn] = 0;
867  //use end of L for temporary space
868  int *indexL = &indexRowL[lSave];
869  CoinFactorizationDouble *multipliersL = &elementL[lSave];
870 
871  //adjust
872  int j;
873 
874  for ( j = 0; j < numberInPivotColumn; j++ ) {
875  multipliersL[j] *= pivotMultiplier;
876  }
877  //zero out fill
878  CoinBigIndex iErase;
879  for ( iErase = 0; iErase < increment2 * numberInPivotRow;
880  iErase++ ) {
881  workArea2[iErase] = 0;
882  }
883  CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
884  unsigned int *temp2 = workArea2;
885  int * nextColumn = nextColumn_.array();
886 
887  //pack down and move to work
888  int jColumn;
889  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
890  int iColumn = saveColumn[jColumn];
891  CoinBigIndex startColumn = startColumnU[iColumn];
892  CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
893  int iRow = indexRowU[startColumn];
894  CoinFactorizationDouble value = elementU[startColumn];
895  double largest;
896  CoinBigIndex put = startColumn;
897  CoinBigIndex positionLargest = -1;
898  CoinFactorizationDouble thisPivotValue = 0.0;
899 
900  //compress column and find largest not updated
901  bool checkLargest;
902  int mark = markRow[iRow];
903 
904  if ( mark == largeInteger+1 ) {
905  largest = fabs ( value );
906  positionLargest = put;
907  put++;
908  checkLargest = false;
909  } else {
910  //need to find largest
911  largest = 0.0;
912  checkLargest = true;
913  if ( mark != largeInteger ) {
914  //will be updated
915  work[mark] = value;
916  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
917  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
918 
919  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
920  added--;
921  } else {
922  thisPivotValue = value;
923  }
924  }
925  CoinBigIndex i;
926  for ( i = startColumn + 1; i < endColumn; i++ ) {
927  iRow = indexRowU[i];
928  value = elementU[i];
929  int mark = markRow[iRow];
930 
931  if ( mark == largeInteger+1 ) {
932  //keep
933  indexRowU[put] = iRow;
934  elementU[put] = value;
935  if ( checkLargest ) {
936  double absValue = fabs ( value );
937 
938  if ( absValue > largest ) {
939  largest = absValue;
940  positionLargest = put;
941  }
942  }
943  put++;
944  } else if ( mark != largeInteger ) {
945  //will be updated
946  work[mark] = value;
947  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
948  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
949 
950  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
951  added--;
952  } else {
953  thisPivotValue = value;
954  }
955  }
956  //slot in pivot
957  elementU[put] = elementU[startColumn];
958  indexRowU[put] = indexRowU[startColumn];
959  if ( positionLargest == startColumn ) {
960  positionLargest = put; //follow if was largest
961  }
962  put++;
963  elementU[startColumn] = thisPivotValue;
964  indexRowU[startColumn] = pivotRow;
965  //clean up counts
966  startColumn++;
967  numberInColumn[iColumn] = put - startColumn;
968  int * numberInColumnPlus = numberInColumnPlus_.array();
969  numberInColumnPlus[iColumn]++;
970  startColumnU[iColumn]++;
971  //how much space have we got
972  int next = nextColumn[iColumn];
973  CoinBigIndex space;
974 
975  space = startColumnU[next] - put - numberInColumnPlus[next];
976  //assume no zero elements
977  if ( numberInPivotColumn > space ) {
978  //getColumnSpace also moves fixed part
979  if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
980  return false;
981  }
982  //redo starts
983  positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
984  startColumn = startColumnU[iColumn];
985  put = startColumn + numberInColumn[iColumn];
986  }
987  double tolerance = zeroTolerance_;
988 
989  int *nextCount = nextCount_.array();
990  for ( j = 0; j < numberInPivotColumn; j++ ) {
991  value = work[j] - thisPivotValue * multipliersL[j];
992  double absValue = fabs ( value );
993 
994  if ( absValue > tolerance ) {
995  work[j] = 0.0;
996  assert (put<lengthAreaU_);
997  elementU[put] = value;
998  indexRowU[put] = indexL[j];
999  if ( absValue > largest ) {
1000  largest = absValue;
1001  positionLargest = put;
1002  }
1003  put++;
1004  } else {
1005  work[j] = 0.0;
1006  added--;
1007  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1008  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1009 
1010  if ( temp2[word] & ( 1 << bit ) ) {
1011  //take out of row list
1012  iRow = indexL[j];
1013  CoinBigIndex start = startRowU[iRow];
1014  CoinBigIndex end = start + numberInRow[iRow];
1015  CoinBigIndex where = start;
1016 
1017  while ( indexColumnU[where] != iColumn ) {
1018  where++;
1019  } /* endwhile */
1020 #if DEBUG_COIN
1021  if ( where >= end ) {
1022  abort ( );
1023  }
1024 #endif
1025  indexColumnU[where] = indexColumnU[end - 1];
1026  numberInRow[iRow]--;
1027  } else {
1028  //make sure won't be added
1029  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1030  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1031 
1032  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1033  }
1034  }
1035  }
1036  numberInColumn[iColumn] = put - startColumn;
1037  //move largest
1038  if ( positionLargest >= 0 ) {
1039  value = elementU[positionLargest];
1040  iRow = indexRowU[positionLargest];
1041  elementU[positionLargest] = elementU[startColumn];
1042  indexRowU[positionLargest] = indexRowU[startColumn];
1043  elementU[startColumn] = value;
1044  indexRowU[startColumn] = iRow;
1045  }
1046  //linked list for column
1047  if ( nextCount[iColumn + numberRows_] != -2 ) {
1048  //modify linked list
1049  deleteLink ( iColumn + numberRows_ );
1050  addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1051  }
1052  temp2 += increment2;
1053  }
1054  //get space for row list
1055  unsigned int *putBase = workArea2;
1056  int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1057  int i = 0;
1058 
1059  // do linked lists and update counts
1060  while ( bigLoops ) {
1061  bigLoops--;
1062  int bit;
1063  for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
1064  unsigned int *putThis = putBase;
1065  int iRow = indexL[i];
1066 
1067  //get space
1068  int number = 0;
1069  int jColumn;
1070 
1071  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1072  unsigned int test = *putThis;
1073 
1074  putThis += increment2;
1075  test = 1 - ( ( test >> bit ) & 1 );
1076  number += test;
1077  }
1078  int next = nextRow[iRow];
1079  CoinBigIndex space;
1080 
1081  space = startRowU[next] - startRowU[iRow];
1082  number += numberInRow[iRow];
1083  if ( space < number ) {
1084  if ( !getRowSpace ( iRow, number ) ) {
1085  return false;
1086  }
1087  }
1088  // now do
1089  putThis = putBase;
1090  next = nextRow[iRow];
1091  number = numberInRow[iRow];
1092  CoinBigIndex end = startRowU[iRow] + number;
1093  int saveIndex = indexColumnU[startRowU[next]];
1094 
1095  //add in
1096  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1097  unsigned int test = *putThis;
1098 
1099  putThis += increment2;
1100  test = 1 - ( ( test >> bit ) & 1 );
1101  indexColumnU[end] = saveColumn[jColumn];
1102  end += test;
1103  }
1104  //put back next one in case zapped
1105  indexColumnU[startRowU[next]] = saveIndex;
1106  markRow[iRow] = static_cast<T>(largeInteger+1);
1107  number = end - startRowU[iRow];
1108  numberInRow[iRow] = number;
1109  deleteLink ( iRow );
1110  addLink ( iRow, number );
1111  }
1112  putBase++;
1113  } /* endwhile */
1114  int bit;
1115 
1116  for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
1117  unsigned int *putThis = putBase;
1118  int iRow = indexL[i];
1119 
1120  //get space
1121  int number = 0;
1122  int jColumn;
1123 
1124  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1125  unsigned int test = *putThis;
1126 
1127  putThis += increment2;
1128  test = 1 - ( ( test >> bit ) & 1 );
1129  number += test;
1130  }
1131  int next = nextRow[iRow];
1132  CoinBigIndex space;
1133 
1134  space = startRowU[next] - startRowU[iRow];
1135  number += numberInRow[iRow];
1136  if ( space < number ) {
1137  if ( !getRowSpace ( iRow, number ) ) {
1138  return false;
1139  }
1140  }
1141  // now do
1142  putThis = putBase;
1143  next = nextRow[iRow];
1144  number = numberInRow[iRow];
1145  CoinBigIndex end = startRowU[iRow] + number;
1146  int saveIndex;
1147 
1148  saveIndex = indexColumnU[startRowU[next]];
1149 
1150  //add in
1151  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1152  unsigned int test = *putThis;
1153 
1154  putThis += increment2;
1155  test = 1 - ( ( test >> bit ) & 1 );
1156 
1157  indexColumnU[end] = saveColumn[jColumn];
1158  end += test;
1159  }
1160  indexColumnU[startRowU[next]] = saveIndex;
1161  markRow[iRow] = static_cast<T>(largeInteger+1);
1162  number = end - startRowU[iRow];
1163  numberInRow[iRow] = number;
1164  deleteLink ( iRow );
1165  addLink ( iRow, number );
1166  }
1167  markRow[pivotRow] = static_cast<T>(largeInteger+1);
1168  //modify linked list for pivots
1169  deleteLink ( pivotRow );
1170  deleteLink ( pivotColumn + numberRows_ );
1171  totalElements_ += added;
1172  return true;
1173 }
1174 
1175  /********************************* END LARGE TEMPLATE ********/
1177 protected:
1179 
1182  double pivotTolerance_;
1186 #ifndef COIN_FAST_CODE
1187  double slackValue_;
1189 #else
1190 #ifndef slackValue_
1191 #define slackValue_ -1.0
1192 #endif
1193 #endif
1194  double areaFactor_;
1197  double relaxCheck_;
1232  int status_;
1233 
1238  //int increasingRows_;
1239 
1244 
1247 
1250 
1253 
1257 
1260 
1263 
1266 
1269 
1272 
1275 
1278 
1281 
1284 
1287 
1290 
1293 
1296 
1299 
1302 
1305 
1307  //int baseU_;
1308 
1311 
1314 
1317 
1320 
1323 
1326 
1329 
1332 
1335 
1338 
1341 
1344 
1347 
1350 
1353 
1356 
1359 
1362 
1365 
1368 
1370  double * denseArea_;
1371 
1374 
1377 
1380 
1383 
1386 
1389 
1391  mutable double ftranCountInput_;
1392  mutable double ftranCountAfterL_;
1393  mutable double ftranCountAfterR_;
1394  mutable double ftranCountAfterU_;
1395  mutable double btranCountInput_;
1396  mutable double btranCountAfterU_;
1397  mutable double btranCountAfterR_;
1398  mutable double btranCountAfterL_;
1399 
1401  mutable int numberFtranCounts_;
1402  mutable int numberBtranCounts_;
1403 
1411 
1413  mutable bool collectStatistics_;
1414 
1417 
1420 
1423 
1426 
1429 
1435  int biasLU_;
1443 };
1444 // Dense coding
1445 #ifdef COIN_HAS_LAPACK
1446 #define DENSE_CODE 1
1447 /* Type of Fortran integer translated into C */
1448 #ifndef ipfint
1449 //typedef ipfint FORTRAN_INTEGER_TYPE ;
1450 typedef int ipfint;
1451 typedef const int cipfint;
1452 #endif
1453 #endif
1454 #endif
1455 // Extra for ugly include
1456 #ifdef UGLY_COIN_FACTOR_CODING
1457 #define FAC_UNSET (FAC_SET+1)
1458 {
1459  goodPivot=false;
1460  //store pivot columns (so can easily compress)
1461  CoinBigIndex startColumnThis = startColumn[iPivotColumn];
1462  CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
1463  int put = 0;
1464  CoinBigIndex startRowThis = startRow[iPivotRow];
1465  CoinBigIndex endRow = startRowThis + numberDoRow + 1;
1466  if ( pivotColumnPosition < 0 ) {
1467  for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1468  int iColumn = indexColumn[pivotColumnPosition];
1469  if ( iColumn != iPivotColumn ) {
1470  saveColumn[put++] = iColumn;
1471  } else {
1472  break;
1473  }
1474  }
1475  } else {
1476  for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
1477  saveColumn[put++] = indexColumn[i];
1478  }
1479  }
1480  assert (pivotColumnPosition<endRow);
1481  assert (indexColumn[pivotColumnPosition]==iPivotColumn);
1482  pivotColumnPosition++;
1483  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1484  saveColumn[put++] = indexColumn[pivotColumnPosition];
1485  }
1486  //take out this bit of indexColumn
1487  int next = nextRow[iPivotRow];
1488  int last = lastRow[iPivotRow];
1489 
1490  nextRow[last] = next;
1491  lastRow[next] = last;
1492  nextRow[iPivotRow] = numberGoodU_; //use for permute
1493  lastRow[iPivotRow] = -2;
1494  numberInRow[iPivotRow] = 0;
1495  //store column in L, compress in U and take column out
1496  CoinBigIndex l = lengthL_;
1497  // **** HORRID coding coming up but a goto seems best!
1498  {
1499  if ( l + numberDoColumn > lengthAreaL_ ) {
1500  //need more memory
1501  if ((messageLevel_&4)!=0)
1502  printf("more memory needed in middle of invert\n");
1503  goto BAD_PIVOT;
1504  }
1505  //l+=currentAreaL_->elementByColumn-elementL;
1506  CoinBigIndex lSave = l;
1507 
1508  CoinBigIndex * startColumnL = startColumnL_.array();
1509  startColumnL[numberGoodL_] = l; //for luck and first time
1510  numberGoodL_++;
1511  startColumnL[numberGoodL_] = l + numberDoColumn;
1512  lengthL_ += numberDoColumn;
1513  if ( pivotRowPosition < 0 ) {
1514  for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1515  int iRow = indexRow[pivotRowPosition];
1516  if ( iRow != iPivotRow ) {
1517  indexRowL[l] = iRow;
1518  elementL[l] = element[pivotRowPosition];
1519  markRow[iRow] = l - lSave;
1520  l++;
1521  //take out of row list
1522  CoinBigIndex start = startRow[iRow];
1523  CoinBigIndex end = start + numberInRow[iRow];
1524  CoinBigIndex where = start;
1525 
1526  while ( indexColumn[where] != iPivotColumn ) {
1527  where++;
1528  } /* endwhile */
1529 #if DEBUG_COIN
1530  if ( where >= end ) {
1531  abort ( );
1532  }
1533 #endif
1534  indexColumn[where] = indexColumn[end - 1];
1535  numberInRow[iRow]--;
1536  } else {
1537  break;
1538  }
1539  }
1540  } else {
1541  CoinBigIndex i;
1542 
1543  for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
1544  int iRow = indexRow[i];
1545 
1546  markRow[iRow] = l - lSave;
1547  indexRowL[l] = iRow;
1548  elementL[l] = element[i];
1549  l++;
1550  //take out of row list
1551  CoinBigIndex start = startRow[iRow];
1552  CoinBigIndex end = start + numberInRow[iRow];
1553  CoinBigIndex where = start;
1554 
1555  while ( indexColumn[where] != iPivotColumn ) {
1556  where++;
1557  } /* endwhile */
1558 #if DEBUG_COIN
1559  if ( where >= end ) {
1560  abort ( );
1561  }
1562 #endif
1563  indexColumn[where] = indexColumn[end - 1];
1564  numberInRow[iRow]--;
1565  assert (numberInRow[iRow]>=0);
1566  }
1567  }
1568  assert (pivotRowPosition<endColumn);
1569  assert (indexRow[pivotRowPosition]==iPivotRow);
1570  CoinFactorizationDouble pivotElement = element[pivotRowPosition];
1571  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1572 
1573  pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1574  pivotRowPosition++;
1575  for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1576  int iRow = indexRow[pivotRowPosition];
1577 
1578  markRow[iRow] = l - lSave;
1579  indexRowL[l] = iRow;
1580  elementL[l] = element[pivotRowPosition];
1581  l++;
1582  //take out of row list
1583  CoinBigIndex start = startRow[iRow];
1584  CoinBigIndex end = start + numberInRow[iRow];
1585  CoinBigIndex where = start;
1586 
1587  while ( indexColumn[where] != iPivotColumn ) {
1588  where++;
1589  } /* endwhile */
1590 #if DEBUG_COIN
1591  if ( where >= end ) {
1592  abort ( );
1593  }
1594 #endif
1595  indexColumn[where] = indexColumn[end - 1];
1596  numberInRow[iRow]--;
1597  assert (numberInRow[iRow]>=0);
1598  }
1599  markRow[iPivotRow] = FAC_SET;
1600  //compress pivot column (move pivot to front including saved)
1601  numberInColumn[iPivotColumn] = 0;
1602  //use end of L for temporary space
1603  int *indexL = &indexRowL[lSave];
1604  CoinFactorizationDouble *multipliersL = &elementL[lSave];
1605 
1606  //adjust
1607  int j;
1608 
1609  for ( j = 0; j < numberDoColumn; j++ ) {
1610  multipliersL[j] *= pivotMultiplier;
1611  }
1612  //zero out fill
1613  CoinBigIndex iErase;
1614  for ( iErase = 0; iErase < increment2 * numberDoRow;
1615  iErase++ ) {
1616  workArea2[iErase] = 0;
1617  }
1618  CoinBigIndex added = numberDoRow * numberDoColumn;
1619  unsigned int *temp2 = workArea2;
1620  int * nextColumn = nextColumn_.array();
1621 
1622  //pack down and move to work
1623  int jColumn;
1624  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1625  int iColumn = saveColumn[jColumn];
1626  CoinBigIndex startColumnThis = startColumn[iColumn];
1627  CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
1628  int iRow = indexRow[startColumnThis];
1629  CoinFactorizationDouble value = element[startColumnThis];
1630  double largest;
1631  CoinBigIndex put = startColumnThis;
1632  CoinBigIndex positionLargest = -1;
1633  CoinFactorizationDouble thisPivotValue = 0.0;
1634 
1635  //compress column and find largest not updated
1636  bool checkLargest;
1637  int mark = markRow[iRow];
1638 
1639  if ( mark == FAC_UNSET ) {
1640  largest = fabs ( value );
1641  positionLargest = put;
1642  put++;
1643  checkLargest = false;
1644  } else {
1645  //need to find largest
1646  largest = 0.0;
1647  checkLargest = true;
1648  if ( mark != FAC_SET ) {
1649  //will be updated
1650  workArea[mark] = value;
1651  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1652  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1653 
1654  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1655  added--;
1656  } else {
1657  thisPivotValue = value;
1658  }
1659  }
1660  CoinBigIndex i;
1661  for ( i = startColumnThis + 1; i < endColumn; i++ ) {
1662  iRow = indexRow[i];
1663  value = element[i];
1664  int mark = markRow[iRow];
1665 
1666  if ( mark == FAC_UNSET ) {
1667  //keep
1668  indexRow[put] = iRow;
1669  element[put] = value;
1670  if ( checkLargest ) {
1671  double absValue = fabs ( value );
1672 
1673  if ( absValue > largest ) {
1674  largest = absValue;
1675  positionLargest = put;
1676  }
1677  }
1678  put++;
1679  } else if ( mark != FAC_SET ) {
1680  //will be updated
1681  workArea[mark] = value;
1682  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1683  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1684 
1685  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1686  added--;
1687  } else {
1688  thisPivotValue = value;
1689  }
1690  }
1691  //slot in pivot
1692  element[put] = element[startColumnThis];
1693  indexRow[put] = indexRow[startColumnThis];
1694  if ( positionLargest == startColumnThis ) {
1695  positionLargest = put; //follow if was largest
1696  }
1697  put++;
1698  element[startColumnThis] = thisPivotValue;
1699  indexRow[startColumnThis] = iPivotRow;
1700  //clean up counts
1701  startColumnThis++;
1702  numberInColumn[iColumn] = put - startColumnThis;
1703  int * numberInColumnPlus = numberInColumnPlus_.array();
1704  numberInColumnPlus[iColumn]++;
1705  startColumn[iColumn]++;
1706  //how much space have we got
1707  int next = nextColumn[iColumn];
1708  CoinBigIndex space;
1709 
1710  space = startColumn[next] - put - numberInColumnPlus[next];
1711  //assume no zero elements
1712  if ( numberDoColumn > space ) {
1713  //getColumnSpace also moves fixed part
1714  if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
1715  goto BAD_PIVOT;
1716  }
1717  //redo starts
1718  positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
1719  startColumnThis = startColumn[iColumn];
1720  put = startColumnThis + numberInColumn[iColumn];
1721  }
1722  double tolerance = zeroTolerance_;
1723 
1724  int *nextCount = nextCount_.array();
1725  for ( j = 0; j < numberDoColumn; j++ ) {
1726  value = workArea[j] - thisPivotValue * multipliersL[j];
1727  double absValue = fabs ( value );
1728 
1729  if ( absValue > tolerance ) {
1730  workArea[j] = 0.0;
1731  element[put] = value;
1732  indexRow[put] = indexL[j];
1733  if ( absValue > largest ) {
1734  largest = absValue;
1735  positionLargest = put;
1736  }
1737  put++;
1738  } else {
1739  workArea[j] = 0.0;
1740  added--;
1741  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1742  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1743 
1744  if ( temp2[word] & ( 1 << bit ) ) {
1745  //take out of row list
1746  iRow = indexL[j];
1747  CoinBigIndex start = startRow[iRow];
1748  CoinBigIndex end = start + numberInRow[iRow];
1749  CoinBigIndex where = start;
1750 
1751  while ( indexColumn[where] != iColumn ) {
1752  where++;
1753  } /* endwhile */
1754 #if DEBUG_COIN
1755  if ( where >= end ) {
1756  abort ( );
1757  }
1758 #endif
1759  indexColumn[where] = indexColumn[end - 1];
1760  numberInRow[iRow]--;
1761  } else {
1762  //make sure won't be added
1763  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1764  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1765 
1766  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1767  }
1768  }
1769  }
1770  numberInColumn[iColumn] = put - startColumnThis;
1771  //move largest
1772  if ( positionLargest >= 0 ) {
1773  value = element[positionLargest];
1774  iRow = indexRow[positionLargest];
1775  element[positionLargest] = element[startColumnThis];
1776  indexRow[positionLargest] = indexRow[startColumnThis];
1777  element[startColumnThis] = value;
1778  indexRow[startColumnThis] = iRow;
1779  }
1780  //linked list for column
1781  if ( nextCount[iColumn + numberRows_] != -2 ) {
1782  //modify linked list
1783  deleteLink ( iColumn + numberRows_ );
1784  addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1785  }
1786  temp2 += increment2;
1787  }
1788  //get space for row list
1789  unsigned int *putBase = workArea2;
1790  int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1791  int i = 0;
1792 
1793  // do linked lists and update counts
1794  while ( bigLoops ) {
1795  bigLoops--;
1796  int bit;
1797  for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
1798  unsigned int *putThis = putBase;
1799  int iRow = indexL[i];
1800 
1801  //get space
1802  int number = 0;
1803  int jColumn;
1804 
1805  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1806  unsigned int test = *putThis;
1807 
1808  putThis += increment2;
1809  test = 1 - ( ( test >> bit ) & 1 );
1810  number += test;
1811  }
1812  int next = nextRow[iRow];
1813  CoinBigIndex space;
1814 
1815  space = startRow[next] - startRow[iRow];
1816  number += numberInRow[iRow];
1817  if ( space < number ) {
1818  if ( !getRowSpace ( iRow, number ) ) {
1819  goto BAD_PIVOT;
1820  }
1821  }
1822  // now do
1823  putThis = putBase;
1824  next = nextRow[iRow];
1825  number = numberInRow[iRow];
1826  CoinBigIndex end = startRow[iRow] + number;
1827  int saveIndex = indexColumn[startRow[next]];
1828 
1829  //add in
1830  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1831  unsigned int test = *putThis;
1832 
1833  putThis += increment2;
1834  test = 1 - ( ( test >> bit ) & 1 );
1835  indexColumn[end] = saveColumn[jColumn];
1836  end += test;
1837  }
1838  //put back next one in case zapped
1839  indexColumn[startRow[next]] = saveIndex;
1840  markRow[iRow] = FAC_UNSET;
1841  number = end - startRow[iRow];
1842  numberInRow[iRow] = number;
1843  deleteLink ( iRow );
1844  addLink ( iRow, number );
1845  }
1846  putBase++;
1847  } /* endwhile */
1848  int bit;
1849 
1850  for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
1851  unsigned int *putThis = putBase;
1852  int iRow = indexL[i];
1853 
1854  //get space
1855  int number = 0;
1856  int jColumn;
1857 
1858  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1859  unsigned int test = *putThis;
1860 
1861  putThis += increment2;
1862  test = 1 - ( ( test >> bit ) & 1 );
1863  number += test;
1864  }
1865  int next = nextRow[iRow];
1866  CoinBigIndex space;
1867 
1868  space = startRow[next] - startRow[iRow];
1869  number += numberInRow[iRow];
1870  if ( space < number ) {
1871  if ( !getRowSpace ( iRow, number ) ) {
1872  goto BAD_PIVOT;
1873  }
1874  }
1875  // now do
1876  putThis = putBase;
1877  next = nextRow[iRow];
1878  number = numberInRow[iRow];
1879  CoinBigIndex end = startRow[iRow] + number;
1880  int saveIndex;
1881 
1882  saveIndex = indexColumn[startRow[next]];
1883 
1884  //add in
1885  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1886  unsigned int test = *putThis;
1887 
1888  putThis += increment2;
1889  test = 1 - ( ( test >> bit ) & 1 );
1890 
1891  indexColumn[end] = saveColumn[jColumn];
1892  end += test;
1893  }
1894  indexColumn[startRow[next]] = saveIndex;
1895  markRow[iRow] = FAC_UNSET;
1896  number = end - startRow[iRow];
1897  numberInRow[iRow] = number;
1898  deleteLink ( iRow );
1899  addLink ( iRow, number );
1900  }
1901  markRow[iPivotRow] = FAC_UNSET;
1902  //modify linked list for pivots
1903  deleteLink ( iPivotRow );
1904  deleteLink ( iPivotColumn + numberRows_ );
1905  totalElements_ += added;
1906  goodPivot= true;
1907  // **** UGLY UGLY UGLY
1908  }
1909  BAD_PIVOT:
1910 
1911  ;
1912 }
1913 #undef FAC_UNSET
1914 #endif
int CoinBigIndex
bool getColumnSpace(int iColumn, int extraNeeded)
Gets space for one Column with given length, may have to do compression (returns True if successful)...
int numberGoodU_
Number factorized in U (not row singletons)
double slackValue() const
Whether slack value is +1 or -1.
bool collectStatistics_
For statistics.
CoinBigIndex * startColumnL() const
Start of each column in L.
void replaceColumnU(CoinIndexedVector *regionSparse, CoinBigIndex *deleted, int internalPivotRow)
Combines BtranU and delete elements If deleted is NULL then delete elements otherwise store where ele...
int addColumn(CoinBigIndex numberElements, int indicesRow[], double elements[])
Adds one Column to basis, can increase size of basis.
int factorDense()
Does dense phase of factorization return code is <0 error, 0= finished.
double ftranCountAfterL_
Pivot tolerance.
CoinBigIndex totalElements_
Number of elements in U (to go) or while iterating total overall.
CoinBigIndex lengthAreaU() const
Returns length of U area.
int maximumColumnsExtra()
Maximum number of Columns after iterating.
CoinBigIndex numberElementsU() const
Returns number in U area.
CoinBigIndexArrayWithLength startRowL_
Start of each row in L.
void updateColumnTransposeLSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparse (by Row)
int numberBtranCounts_
Pivot tolerance.
void updateColumnTransposeUSparsish(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) when sparsish, assumes index is sorted i.e.
CoinBigIndex getColumnSpaceIterate(int iColumn, double value, int iRow)
getColumnSpaceIterate.
int addRow(CoinBigIndex numberElements, int indicesColumn[], double elements[])
Adds one Row to basis, can increase size of basis.
void updateColumnTransposePFI(CoinIndexedVector *region) const
Updates part of column transpose PFI (BTRAN) (before rest)
void setBiasLU(int value)
Returns number of dense rows.
int numberForrestTomlin() const
Length of FT vector.
void gutsOfCopy(const CoinFactorization &other)
See if worth going sparse.
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 * indexColumnL() const
Index of column in row for L.
double zeroTolerance_
Zero tolerance.
int sparseThreshold_
Below this use sparse technology - if 0 then no L row copy.
CoinIntArrayWithLength saveColumn_
Columns left to do in a single pivot.
int add(CoinBigIndex numberElements, int indicesRow[], int indicesColumn[], double elements[])
Adds given elements to Basis and updates factorization, can increase size of basis.
void clearArrays()
Get rid of all memory.
CoinUnsignedIntArrayWithLength workArea2_
Second work area.
void updateColumnTransposeUDensish(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) when densish, assumes index is sorted i.e.
int * pivotColumn() const
Returns address of pivotColumn region (also used for permuting)
double ftranCountAfterU_
Pivot tolerance.
void checkConsistency()
Checks that row and column copies look OK.
CoinBigIndex lengthR_
Length of R stuff.
bool pivot(int pivotRow, int pivotColumn, CoinBigIndex pivotRowPosition, CoinBigIndex pivotColumnPosition, CoinFactorizationDouble work[], unsigned int workArea2[], int increment2, T markRow[], int largeInteger)
Gets space for a factorization, called by constructors.
int numberTrials_
0 - no increasing rows - no permutations, 1 - no increasing rows but permutations 2 - increasing rows...
CoinBigIndex baseL() const
Base of L.
int deleteColumn(int Row)
Deletes one Column from basis, returns rank.
CoinBigIndex lengthAreaU_
Length of area reserved for U.
void preProcess(int state, int possibleDuplicates=-1)
PreProcesses raw triplet data.
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 addLink(int index, int count)
Adds a link in chain of equal counts.
CoinBigIndex numberL() const
Number in L.
CoinIntArrayWithLength indexColumnU_
Base address for U (may change)
CoinFactorizationDouble * elementR_
Elements of R.
CoinIntArrayWithLength nextRow_
Next Row in memory order.
bool pivotRowSingleton(int pivotRow, int pivotColumn)
Does one pivot on Row Singleton in factorization.
double btranCountAfterR_
Pivot tolerance.
int numberColumns() const
Total number of columns in factorization.
void gutsOfDestructor(int type=1)
The real work of constructors etc 0 just scalars, 1 bit normal.
#define COINFACTORIZATION_BITS_PER_INT
int numberColumns_
Number of Columns in factorization.
int * indexRowU() const
Row indices of U.
int biggerDimension_
Larger of row and column size.
double CoinFactorizationDouble
Definition: CoinFinite.hpp:86
int sparseThreshold2_
And one for "sparsish".
CoinBigIndex lengthU_
Base of U is always 0.
bool getRowSpace(int iRow, int extraNeeded)
Gets space for one Row with given length, may have to do compression (returns True if successful)...
void updateColumnTransposeRSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR) when sparse.
CoinBigIndex numberL_
Number in L.
double ftranCountAfterR_
Pivot tolerance.
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 indexColumnL_
Index of column in row for L.
int sparseThreshold() const
get sparse threshold
int deleteRow(int Row)
Deletes one Row from basis, returns rank.
int * array() const
Get Array.
int updateColumnUDensish(double *COIN_RESTRICT region, int *COIN_RESTRICT regionIndex) const
Updates part of column (FTRANU)
int * indexRowL() const
Row indices of L.
double ftranAverageAfterU_
Pivot tolerance.
void setPersistenceFlag(int value)
Returns number of dense rows.
CoinFactorizationDoubleArrayWithLength elementU_
Elements of U.
CoinIntArrayWithLength sparse_
Sparse regions.
CoinBigIndex numberElementsL() const
Returns number in L area.
void updateColumnTransposeUByColumn(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU) by column assumes index is sorted i.e.
CoinFactorizationDoubleArrayWithLength workArea_
First work area.
Sparse Matrix Base Class.
CoinBigIndex maximumU_
Maximum space used in U.
CoinIntArrayWithLength permute_
Permutation vector for pivot row order.
CoinBigIndex numberCompressions_
Number of compressions done.
double ftranAverageAfterL_
While these are average ratios collected over last period.
bool forrestTomlin() const
true if Forrest Tomlin update, false if PFI
void updateColumnUSparse(CoinIndexedVector *regionSparse, int *indexIn) const
Updates part of column (FTRANU) when sparse.
void setNumberRows(int value)
Set number of Rows after factorization.
CoinIntArrayWithLength lastCount_
Previous Row/Column with count.
CoinBigIndex numberElementsR() const
Returns number in R area.
bool reorderU()
Reorders U so contiguous and in order (if there is space) Returns true if it could.
CoinFactorization()
Default constructor.
bool collectStatistics() const
For statistics.
void updateColumnTransposeUSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANU) when sparse, assumes index is sorted i.e.
double adjustedAreaFactor() const
Returns areaFactor but adjusted for dense.
bool spaceForForrestTomlin() const
True if FT update and space.
CoinIntArrayWithLength numberInColumnPlus_
Number in each Column including pivoted.
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...
void areaFactor(double value)
Returns status.
CoinBigIndexArrayWithLength startRowU_
Start of each Row as pointer.
bool pivotOneOtherRow(int pivotRow, int pivotColumn)
Pivots when just one other row so faster?
void goSparse()
makes a row copy of L for speed and to allow very sparse problems
int numberColumnsExtra_
Number of Columns after iterating.
double ftranCountInput_
Below are all to collect.
CoinBigIndex factorElements_
Number of elements after factorization.
CoinBigIndex lengthAreaL_
Length of area reserved for L.
void updateColumnUSparsish(CoinIndexedVector *regionSparse, int *indexIn) const
Updates part of column (FTRANU) when sparsish.
CoinFactorizationDouble * elementByRowL() const
Elements in L (row copy)
void updateColumnLDensish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when densish.
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 maximumPivots_
Maximum number of pivots before factorization.
int maximumColumnsExtra_
Maximum number of Columns after iterating.
CoinFactorizationDoubleArrayWithLength pivotRegion_
Inverses of pivot values.
CoinIntArrayWithLength pivotColumn_
Pivot order for each Column.
void updateColumnTransposeU(CoinIndexedVector *region, int smallestIndex) const
Updates part of column transpose (BTRANU), assumes index is sorted i.e.
CoinBigIndexArrayWithLength startColumnR_
Start of columns for R.
CoinFactorizationDouble * elementU() const
Elements of U.
void updateColumnTransposeLDensish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when densish by column.
CoinIntArrayWithLength indexRowL_
Row indices of L.
double btranAverageAfterR_
Pivot tolerance.
double * denseArea_
Dense area.
CoinIntArrayWithLength nextColumn_
Next Column in memory order.
void gutsOfInitialize(int type)
1 bit - tolerances etc, 2 more, 4 dummy arrays
int maximumPivots() const
Maximum number of pivots between factorizations.
int biasLU() const
L to U bias 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias.
bool getColumnSpaceIterateR(int iColumn, double value, int iRow)
getColumnSpaceIterateR.
void updateColumnTransposeL(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL)
CoinIntArrayWithLength numberInRow_
Number in each Row.
void setNumberElementsU(CoinBigIndex value)
Setss number in U area.
void sort() const
Debug - sort so can compare.
CoinFactorizationDouble * array() const
Get Array.
#define COINFACTORIZATION_MASK_PER_INT
void updateColumnTransposeR(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR)
CoinBigIndex * array() const
Get Array.
CoinFactorizationDoubleArrayWithLength elementL_
Elements of L.
friend void CoinFactorizationUnitTest(const std::string &mpsDir)
#define COINFACTORIZATION_SHIFT_PER_INT
int updateColumnFT(CoinIndexedVector *regionSparse, CoinIndexedVector *regionSparse2)
Updates one column (FTRAN) from regionSparse2 Tries to do FT update number returned is negative if no...
int * indexRowR_
Row indices for R.
int maximumRowsExtra() const
Maximum of Rows after iterating.
CoinBigIndex lengthAreaR_
length of area reserved for R
double slackValue_
Whether slack value is +1 or -1.
double areaFactor() const
Whether larger areas needed.
double btranCountAfterL_
Pivot tolerance.
void checkSparse()
See if worth going sparse.
int * numberInColumn() const
Number of entries in each column.
int * permute() const
Returns address of permute region.
int numberSlacks_
Number of slacks at beginning of U.
void resetStatistics()
Reset all sparsity etc statistics.
Indexed Vector.
int * numberInRow() const
Number of entries in each row.
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.
CoinBigIndexArrayWithLength startColumnU_
Start of each column in U.
void setCollectStatistics(bool onOff) const
For statistics.
void updateColumnRFT(CoinIndexedVector *region, int *indexIn)
Updates part of column (FTRANR) with FT update.
CoinIntArrayWithLength firstCount_
First Row/Column with count of k, can tell which by offset - Rows then Columns.
void getAreas(int numberRows, int numberColumns, CoinBigIndex maximumL, CoinBigIndex maximumU)
Gets space for a factorization, called by constructors.
int numberElements() const
Total number of elements in factorization.
double getAccuracyCheck() const
Returns status.
double btranCountInput_
Pivot tolerance.
#define COIN_RESTRICT
Definition: CoinFinite.hpp:45
double btranCountAfterU_
Pivot tolerance.
int numberU_
Number in U.
int factorSparseSmall()
Does sparse phase of factorization (for smaller problems) return code is <0 error, 0= finished.
CoinIntArrayWithLength markRow_
Marks rows to be updated.
~CoinFactorization()
Destructor.
int factorize(const CoinPackedMatrix &matrix, int rowIsBasic[], int columnIsBasic[], double areaFactor=0.0)
When part of LP - given by basic variables.
CoinFactorizationDouble * pivotRegion() const
Returns address of pivot region.
void cleanup()
Cleans up at end of factorization.
CoinIntArrayWithLength pivotRowL_
Pivots for L.
CoinIntArrayWithLength lastRow_
Previous Row in memory order.
int numberPivots_
Number pivots since last factorization.
int * permuteBack() const
Returns address of permuteBack region.
CoinBigIndex baseL_
Base of L.
int denseThreshold_
Dense threshold.
CoinBigIndexArrayWithLength startColumnL_
Start of each column in L.
int numberDense() const
Returns number of dense rows.
CoinIntArrayWithLength nextCount_
Next Row/Column with count.
CoinBigIndex * startColumnU() const
Start of each column in U.
int factor()
Does most of factorization.
CoinBigIndex * startRowL() const
Start of each row in L.
void separateLinks(int count, bool rowsFirst)
Separate out links with same row/column count.
void updateColumnL(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL)
int status_
Status of factorization.
void updateColumnTransposeLByRow(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when densish by row.
int numberRowsExtra() const
Number of Rows after iterating.
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 factorSparseLarge()
Does sparse phase of factorization (for larger problems) return code is <0 error, 0= finished...
int numberFtranCounts_
We can roll over factorizations.
int pivots() const
Returns number of pivots since factorization.
void setPivots(int value)
Sets number of pivots since factorization.
CoinIntArrayWithLength numberInColumn_
Number in each Column.
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...
void updateColumnR(CoinIndexedVector *region) const
Updates part of column (FTRANR) without FT update.
double maximumCoefficient() const
Returns maximum absolute value in factorization.
double areaFactor_
How much to multiply areas by.
void updateColumnU(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANU)
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 numberRows_
Number of Rows in factorization.
int biasLU_
L to U bias 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias.
double btranAverageAfterU_
Pivot tolerance.
int factorSparse()
Does sparse phase of factorization return code is <0 error, 0= finished.
int status() const
Returns status.
CoinBigIndex numberCompressions() const
Number of compressions done.
int numberGoodColumns() const
Number of good columns in factorization.
void setDenseThreshold(int value)
Sets dense threshold.
CoinBigIndexArrayWithLength convertRowToColumnU_
Converts rows to columns in U.
int checkPivot(double saveFromU, double oldPivot) const
Returns accuracy status of replaceColumn returns 0=OK, 1=Probably OK, 2=singular. ...
CoinBigIndex * version.
void updateColumnPFI(CoinIndexedVector *regionSparse) const
Updates part of column PFI (FTRAN) (after rest)
int numberRows() const
Number of Rows after factorization.
void emptyRows(int numberToEmpty, const int which[])
Takes out all entries for given rows.
This deals with Factorization and Updates.
int numberGoodL_
Number factorized in L.
double pivotTolerance() const
Pivot tolerance.
double btranAverageAfterL_
Pivot tolerance.
CoinIntArrayWithLength permuteBack_
DePermutation vector for pivot row order.
int restoreFactorization(const char *file, bool factor=false)
Debug - restore from file - 0 if no error on file.
double pivotTolerance_
Pivot tolerance.
int factorizePart1(int numberRows, int numberColumns, CoinBigIndex 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 updateColumnTransposeRDensish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANR) when dense.
double relaxCheck_
Relax check on accuracy in replaceColumn.
int denseThreshold() const
Gets dense threshold.
CoinBigIndex lengthAreaL() const
Returns length of L area.
void almostDestructor()
Delete all stuff (leaves as after CoinFactorization())
CoinFactorizationDouble * version.
void setStatus(int value)
Sets status.
CoinIntArrayWithLength indexRowU_
Row indices of U.
bool doForrestTomlin_
true if Forrest Tomlin update, false if PFI
CoinIntArrayWithLength lastColumn_
Previous Column in memory order.
void relaxAccuracyCheck(double value)
Allows change of pivot accuracy check 1.0 == none >1.0 relaxed.
CoinBigIndex lengthL_
Length of L.
int messageLevel_
Detail in messages.
int saveFactorization(const char *file) const
Debug - save on file - 0 if no error.
void show_self() const
Debug show object (shows one representation)
int maximumRowsExtra_
Maximum number of Rows after iterating.
double zeroTolerance() const
Zero tolerance.
int numberDense_
Number of dense rows.
int messageLevel() const
Level of detail of messages.
int factorizePart2(int permutation[], int exactNumberElements)
This is part two of factorization Arrays belong to factorization and were returned by part 1 If statu...
double ftranAverageAfterR_
Pivot tolerance.
int numberRowsExtra_
Number of Rows after iterating.
void updateColumnLSparse(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when sparse.
CoinIntArrayWithLength pivotColumnBack_
Inverse Pivot order for each Column.
void deleteLink(int index)
Deletes a link in chain of equal counts.
CoinFactorizationDoubleArrayWithLength elementByRowL_
Elements in L (row copy)
bool pivotColumnSingleton(int pivotRow, int pivotColumn)
Does one pivot on Column Singleton in factorization.
int numberR_
Number in R.
int * densePermute_
Dense permutation.
CoinFactorization & operator=(const CoinFactorization &other)
= copy
void updateColumnLSparsish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when sparsish.
int * pivotColumnBack() const
Returns address of pivotColumnBack region (also used for permuting) Now uses firstCount to save memor...
double conditionNumber() const
Condition number - product of pivots after factorization.
void updateColumnTransposeLSparsish(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparsish 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 setForrestTomlin(bool value)
Returns status.