CoinUtils  2.9.15
 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 1590 2013-04-10 16:48:33Z stefan $ */
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 COIN_ONE_ETA_COPY 100
15 
16 #include <iostream>
17 #include <string>
18 #include <cassert>
19 #include <cstdio>
20 #include <cmath>
21 #include "CoinTypes.hpp"
22 #include "CoinIndexedVector.hpp"
23 
24 class CoinPackedMatrix;
51  friend void CoinFactorizationUnitTest( const std::string & mpsDir );
52 
53 public:
54 
60  CoinFactorization ( const CoinFactorization &other);
61 
65  void almostDestructor();
67  void show_self ( ) const;
69  int saveFactorization (const char * file ) const;
73  int restoreFactorization (const char * file , bool factor=false) ;
75  void sort ( ) const;
79 
89  int factorize ( const CoinPackedMatrix & matrix,
90  int rowIsBasic[], int columnIsBasic[] ,
91  double areaFactor = 0.0 );
102  int factorize ( int numberRows,
103  int numberColumns,
105  CoinBigIndex maximumL,
106  CoinBigIndex maximumU,
107  const int indicesRow[],
108  const int indicesColumn[], const double elements[] ,
109  int permutation[],
110  double areaFactor = 0.0);
115  int factorizePart1 ( int numberRows,
116  int numberColumns,
117  CoinBigIndex estimateNumberElements,
118  int * indicesRow[],
119  int * indicesColumn[],
120  CoinFactorizationDouble * elements[],
121  double areaFactor = 0.0);
128  int factorizePart2 (int permutation[],int exactNumberElements);
130  double conditionNumber() const;
131 
133 
136  inline int status ( ) const {
138  return status_;
139  }
141  inline void setStatus ( int value)
142  { status_=value; }
144  inline int pivots ( ) const {
145  return numberPivots_;
146  }
148  inline void setPivots ( int value )
149  { numberPivots_=value; }
151  inline int *permute ( ) const {
152  return permute_.array();
153  }
155  inline int *pivotColumn ( ) const {
156  return pivotColumn_.array();
157  }
160  return pivotRegion_.array();
161  }
163  inline int *permuteBack ( ) const {
164  return permuteBack_.array();
165  }
168  inline int *pivotColumnBack ( ) const {
169  //return firstCount_.array();
170  return pivotColumnBack_.array();
171  }
173  inline CoinBigIndex * startRowL() const
174  { return startRowL_.array();}
175 
177  inline CoinBigIndex * startColumnL() const
178  { return startColumnL_.array();}
179 
181  inline int * indexColumnL() const
182  { return indexColumnL_.array();}
183 
185  inline int * indexRowL() const
186  { return indexRowL_.array();}
187 
190  { return elementByRowL_.array();}
191 
193  inline int numberRowsExtra ( ) const {
194  return numberRowsExtra_;
195  }
197  inline void setNumberRows(int value)
198  { numberRows_ = value; }
200  inline int numberRows ( ) const {
201  return numberRows_;
202  }
204  inline CoinBigIndex numberL() const
205  { return numberL_;}
206 
208  inline CoinBigIndex baseL() const
209  { return baseL_;}
211  inline int maximumRowsExtra ( ) const {
212  return maximumRowsExtra_;
213  }
215  inline int numberColumns ( ) const {
216  return numberColumns_;
217  }
219  inline int numberElements ( ) const {
220  return totalElements_;
221  }
223  inline int numberForrestTomlin ( ) const {
225  }
227  inline int numberGoodColumns ( ) const {
228  return numberGoodU_;
229  }
231  inline double areaFactor ( ) const {
232  return areaFactor_;
233  }
234  inline void areaFactor ( double value ) {
235  areaFactor_=value;
236  }
238  double adjustedAreaFactor() const;
240  inline void relaxAccuracyCheck(double value)
241  { relaxCheck_ = value;}
242  inline double getAccuracyCheck() const
243  { return relaxCheck_;}
245  inline int messageLevel ( ) const {
246  return messageLevel_ ;
247  }
248  void messageLevel ( int value );
250  inline int maximumPivots ( ) const {
251  return maximumPivots_ ;
252  }
253  void maximumPivots ( int value );
254 
256  inline int denseThreshold() const
257  { return denseThreshold_;}
259  inline void setDenseThreshold(int value)
260  { denseThreshold_ = value;}
262  inline double pivotTolerance ( ) const {
263  return pivotTolerance_ ;
264  }
265  void pivotTolerance ( double value );
267  inline double zeroTolerance ( ) const {
268  return zeroTolerance_ ;
269  }
270  void zeroTolerance ( double value );
271 #ifndef COIN_FAST_CODE
272  inline double slackValue ( ) const {
274  return slackValue_ ;
275  }
276  void slackValue ( double value );
277 #endif
278  double maximumCoefficient() const;
281  inline bool forrestTomlin() const
282  { return doForrestTomlin_;}
283  inline void setForrestTomlin(bool value)
284  { doForrestTomlin_=value;}
286  inline bool spaceForForrestTomlin() const
287  {
289  CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
290  return (space>=0)&&doForrestTomlin_;
291  }
293 
296 
298  inline int numberDense() const
299  { return numberDense_;}
300 
302  inline CoinBigIndex numberElementsU ( ) const {
303  return lengthU_;
304  }
306  inline void setNumberElementsU(CoinBigIndex value)
307  { lengthU_ = value; }
309  inline CoinBigIndex lengthAreaU ( ) const {
310  return lengthAreaU_;
311  }
313  inline CoinBigIndex numberElementsL ( ) const {
314  return lengthL_;
315  }
317  inline CoinBigIndex lengthAreaL ( ) const {
318  return lengthAreaL_;
319  }
321  inline CoinBigIndex numberElementsR ( ) const {
322  return lengthR_;
323  }
326  { return numberCompressions_;}
328  inline int * numberInRow() const
329  { return numberInRow_.array();}
331  inline int * numberInColumn() const
332  { return numberInColumn_.array();}
335  { return elementU_.array();}
337  inline int * indexRowU() const
338  { return indexRowU_.array();}
340  inline CoinBigIndex * startColumnU() const
341  { return startColumnU_.array();}
343  inline int maximumColumnsExtra()
344  { return maximumColumnsExtra_;}
348  inline int biasLU() const
349  { return biasLU_;}
350  inline void setBiasLU(int value)
351  { biasLU_=value;}
357  inline int persistenceFlag() const
358  { return persistenceFlag_;}
359  void setPersistenceFlag(int value);
361 
364 
372  int replaceColumn ( CoinIndexedVector * regionSparse,
373  int pivotRow,
374  double pivotCheck ,
375  bool checkBeforeModifying=false,
376  double acceptablePivot=1.0e-8);
381  void replaceColumnU ( CoinIndexedVector * regionSparse,
382  CoinBigIndex * deleted,
383  int internalPivotRow);
385 
395  int updateColumnFT ( CoinIndexedVector * regionSparse,
396  CoinIndexedVector * regionSparse2);
399  int updateColumn ( CoinIndexedVector * regionSparse,
400  CoinIndexedVector * regionSparse2,
401  bool noPermute=false) const;
407  int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
408  CoinIndexedVector * regionSparse2,
409  CoinIndexedVector * regionSparse3,
410  bool noPermuteRegion3=false) ;
415  int updateColumnTranspose ( CoinIndexedVector * regionSparse,
416  CoinIndexedVector * regionSparse2) const;
418  void goSparse();
420  inline int sparseThreshold ( ) const
421  { return sparseThreshold_;}
423  void sparseThreshold ( int value );
425 
430  inline void clearArrays()
432  { gutsOfDestructor();}
434 
437 
440  int add ( CoinBigIndex numberElements,
441  int indicesRow[],
442  int indicesColumn[], double elements[] );
443 
446  int addColumn ( CoinBigIndex numberElements,
447  int indicesRow[], double elements[] );
448 
451  int addRow ( CoinBigIndex numberElements,
452  int indicesColumn[], double elements[] );
453 
455  int deleteColumn ( int Row );
457  int deleteRow ( int Row );
458 
462  int replaceRow ( int whichRow, int numberElements,
463  const int indicesColumn[], const double elements[] );
465  void emptyRows(int numberToEmpty, const int which[]);
467 
468  void checkSparse();
471  inline bool collectStatistics() const
472  { return collectStatistics_;}
474  inline void setCollectStatistics(bool onOff) const
475  { collectStatistics_ = onOff;}
477  void gutsOfDestructor(int type=1);
479  void gutsOfInitialize(int type);
480  void gutsOfCopy(const CoinFactorization &other);
481 
483  void resetStatistics();
484 
485 
487 
489  void getAreas ( int numberRows,
491  int numberColumns,
492  CoinBigIndex maximumL,
493  CoinBigIndex maximumU );
494 
497  void preProcess ( int state,
498  int possibleDuplicates = -1 );
500  int factor ( );
501 protected:
504  int factorSparse ( );
507  int factorSparseSmall ( );
510  int factorSparseLarge ( );
513  int factorDense ( );
514 
516  bool pivotOneOtherRow ( int pivotRow,
517  int pivotColumn );
519  bool pivotRowSingleton ( int pivotRow,
520  int pivotColumn );
522  bool pivotColumnSingleton ( int pivotRow,
523  int pivotColumn );
524 
529  bool getColumnSpace ( int iColumn,
530  int extraNeeded );
531 
534  bool reorderU();
538  bool getColumnSpaceIterateR ( int iColumn, double value,
539  int iRow);
545  CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
546  int iRow);
550  bool getRowSpace ( int iRow, int extraNeeded );
551 
555  bool getRowSpaceIterate ( int iRow,
556  int extraNeeded );
558  void checkConsistency ( );
560  inline void addLink ( int index, int count ) {
561  int *nextCount = nextCount_.array();
562  int *firstCount = firstCount_.array();
563  int *lastCount = lastCount_.array();
564  int next = firstCount[count];
565  lastCount[index] = -2 - count;
566  if ( next < 0 ) {
567  //first with that count
568  firstCount[count] = index;
569  nextCount[index] = -1;
570  } else {
571  firstCount[count] = index;
572  nextCount[index] = next;
573  lastCount[next] = index;
574  }}
576  inline void deleteLink ( int index ) {
577  int *nextCount = nextCount_.array();
578  int *firstCount = firstCount_.array();
579  int *lastCount = lastCount_.array();
580  int next = nextCount[index];
581  int last = lastCount[index];
582  if ( last >= 0 ) {
583  nextCount[last] = next;
584  } else {
585  int count = -last - 2;
586 
587  firstCount[count] = next;
588  }
589  if ( next >= 0 ) {
590  lastCount[next] = last;
591  }
592  nextCount[index] = -2;
593  lastCount[index] = -2;
594  return;
595  }
597  void separateLinks(int count,bool rowsFirst);
599  void cleanup ( );
600 
602  void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
604  void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
606  void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
608  void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
609 
611  void updateColumnR ( CoinIndexedVector * region ) const;
614  void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
615 
617  void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
618 
620  void updateColumnUSparse ( CoinIndexedVector * regionSparse,
621  int * indexIn) const;
623  void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
624  int * indexIn) const;
626  int updateColumnUDensish ( double * COIN_RESTRICT region,
627  int * COIN_RESTRICT regionIndex) const;
630  int & numberNonZero1,
631  double * COIN_RESTRICT region1,
632  int * COIN_RESTRICT index1,
633  int & numberNonZero2,
634  double * COIN_RESTRICT region2,
635  int * COIN_RESTRICT index2) const;
637  void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
639  void permuteBack ( CoinIndexedVector * regionSparse,
640  CoinIndexedVector * outVector) const;
641 
643  void updateColumnTransposePFI ( CoinIndexedVector * region) const;
647  int smallestIndex) const;
651  int smallestIndex) const;
655  int smallestIndex) const;
658  void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
662  int smallestIndex) const;
663 
665  void updateColumnTransposeR ( CoinIndexedVector * region ) const;
667  void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
669  void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
670 
672  void updateColumnTransposeL ( CoinIndexedVector * region ) const;
674  void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
676  void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
678  void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
680  void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
681 public:
686  int replaceColumnPFI ( CoinIndexedVector * regionSparse,
687  int pivotRow, double alpha);
688 protected:
691  int checkPivot(double saveFromU, double oldPivot) const;
692  /********************************* START LARGE TEMPLATE ********/
693 #ifdef INT_IS_8
694 #define COINFACTORIZATION_BITS_PER_INT 64
695 #define COINFACTORIZATION_SHIFT_PER_INT 6
696 #define COINFACTORIZATION_MASK_PER_INT 0x3f
697 #else
698 #define COINFACTORIZATION_BITS_PER_INT 32
699 #define COINFACTORIZATION_SHIFT_PER_INT 5
700 #define COINFACTORIZATION_MASK_PER_INT 0x1f
701 #endif
702  template <class T> inline bool
703  pivot ( int pivotRow,
704  int pivotColumn,
705  CoinBigIndex pivotRowPosition,
706  CoinBigIndex pivotColumnPosition,
708  unsigned int workArea2[],
709  int increment2,
710  T markRow[] ,
711  int largeInteger)
712 {
713  int *indexColumnU = indexColumnU_.array();
717  int *indexRowU = indexRowU_.array();
718  CoinBigIndex *startRowU = startRowU_.array();
719  int *numberInRow = numberInRow_.array();
721  int *indexRowL = indexRowL_.array();
722  int *saveColumn = saveColumn_.array();
723  int *nextRow = nextRow_.array();
724  int *lastRow = lastRow_.array() ;
725 
726  //store pivot columns (so can easily compress)
727  int numberInPivotRow = numberInRow[pivotRow] - 1;
728  CoinBigIndex startColumn = startColumnU[pivotColumn];
729  int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
730  CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
731  int put = 0;
732  CoinBigIndex startRow = startRowU[pivotRow];
733  CoinBigIndex endRow = startRow + numberInPivotRow + 1;
734 
735  if ( pivotColumnPosition < 0 ) {
736  for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
737  int iColumn = indexColumnU[pivotColumnPosition];
738  if ( iColumn != pivotColumn ) {
739  saveColumn[put++] = iColumn;
740  } else {
741  break;
742  }
743  }
744  } else {
745  for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
746  saveColumn[put++] = indexColumnU[i];
747  }
748  }
749  assert (pivotColumnPosition<endRow);
750  assert (indexColumnU[pivotColumnPosition]==pivotColumn);
751  pivotColumnPosition++;
752  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
753  saveColumn[put++] = indexColumnU[pivotColumnPosition];
754  }
755  //take out this bit of indexColumnU
756  int next = nextRow[pivotRow];
757  int last = lastRow[pivotRow];
758 
759  nextRow[last] = next;
760  lastRow[next] = last;
761  nextRow[pivotRow] = numberGoodU_; //use for permute
762  lastRow[pivotRow] = -2;
763  numberInRow[pivotRow] = 0;
764  //store column in L, compress in U and take column out
766 
767  if ( l + numberInPivotColumn > lengthAreaL_ ) {
768  //need more memory
769  if ((messageLevel_&4)!=0)
770  printf("more memory needed in middle of invert\n");
771  return false;
772  }
773  //l+=currentAreaL_->elementByColumn-elementL;
774  CoinBigIndex lSave = l;
775 
777  startColumnL[numberGoodL_] = l; //for luck and first time
778  numberGoodL_++;
779  startColumnL[numberGoodL_] = l + numberInPivotColumn;
780  lengthL_ += numberInPivotColumn;
781  if ( pivotRowPosition < 0 ) {
782  for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
783  int iRow = indexRowU[pivotRowPosition];
784  if ( iRow != pivotRow ) {
785  indexRowL[l] = iRow;
786  elementL[l] = elementU[pivotRowPosition];
787  markRow[iRow] = static_cast<T>(l - lSave);
788  l++;
789  //take out of row list
790  CoinBigIndex start = startRowU[iRow];
791  CoinBigIndex end = start + numberInRow[iRow];
792  CoinBigIndex where = start;
793 
794  while ( indexColumnU[where] != pivotColumn ) {
795  where++;
796  } /* endwhile */
797 #if DEBUG_COIN
798  if ( where >= end ) {
799  abort ( );
800  }
801 #endif
802  indexColumnU[where] = indexColumnU[end - 1];
803  numberInRow[iRow]--;
804  } else {
805  break;
806  }
807  }
808  } else {
809  CoinBigIndex i;
810 
811  for ( i = startColumn; i < pivotRowPosition; i++ ) {
812  int iRow = indexRowU[i];
813 
814  markRow[iRow] = static_cast<T>(l - lSave);
815  indexRowL[l] = iRow;
816  elementL[l] = elementU[i];
817  l++;
818  //take out of row list
819  CoinBigIndex start = startRowU[iRow];
820  CoinBigIndex end = start + numberInRow[iRow];
821  CoinBigIndex where = start;
822 
823  while ( indexColumnU[where] != pivotColumn ) {
824  where++;
825  } /* endwhile */
826 #if DEBUG_COIN
827  if ( where >= end ) {
828  abort ( );
829  }
830 #endif
831  indexColumnU[where] = indexColumnU[end - 1];
832  numberInRow[iRow]--;
833  assert (numberInRow[iRow]>=0);
834  }
835  }
836  assert (pivotRowPosition<endColumn);
837  assert (indexRowU[pivotRowPosition]==pivotRow);
838  CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
839  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
840 
841  pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
842  pivotRowPosition++;
843  for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
844  int iRow = indexRowU[pivotRowPosition];
845 
846  markRow[iRow] = static_cast<T>(l - lSave);
847  indexRowL[l] = iRow;
848  elementL[l] = elementU[pivotRowPosition];
849  l++;
850  //take out of row list
851  CoinBigIndex start = startRowU[iRow];
852  CoinBigIndex end = start + numberInRow[iRow];
853  CoinBigIndex where = start;
854 
855  while ( indexColumnU[where] != pivotColumn ) {
856  where++;
857  } /* endwhile */
858 #if DEBUG_COIN
859  if ( where >= end ) {
860  abort ( );
861  }
862 #endif
863  indexColumnU[where] = indexColumnU[end - 1];
864  numberInRow[iRow]--;
865  assert (numberInRow[iRow]>=0);
866  }
867  markRow[pivotRow] = static_cast<T>(largeInteger);
868  //compress pivot column (move pivot to front including saved)
869  numberInColumn[pivotColumn] = 0;
870  //use end of L for temporary space
871  int *indexL = &indexRowL[lSave];
872  CoinFactorizationDouble *multipliersL = &elementL[lSave];
873 
874  //adjust
875  int j;
876 
877  for ( j = 0; j < numberInPivotColumn; j++ ) {
878  multipliersL[j] *= pivotMultiplier;
879  }
880  //zero out fill
881  CoinBigIndex iErase;
882  for ( iErase = 0; iErase < increment2 * numberInPivotRow;
883  iErase++ ) {
884  workArea2[iErase] = 0;
885  }
886  CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
887  unsigned int *temp2 = workArea2;
888  int * nextColumn = nextColumn_.array();
889 
890  //pack down and move to work
891  int jColumn;
892  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
893  int iColumn = saveColumn[jColumn];
894  CoinBigIndex startColumn = startColumnU[iColumn];
895  CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
896  int iRow = indexRowU[startColumn];
897  CoinFactorizationDouble value = elementU[startColumn];
898  double largest;
899  CoinBigIndex put = startColumn;
900  CoinBigIndex positionLargest = -1;
901  CoinFactorizationDouble thisPivotValue = 0.0;
902 
903  //compress column and find largest not updated
904  bool checkLargest;
905  int mark = markRow[iRow];
906 
907  if ( mark == largeInteger+1 ) {
908  largest = fabs ( value );
909  positionLargest = put;
910  put++;
911  checkLargest = false;
912  } else {
913  //need to find largest
914  largest = 0.0;
915  checkLargest = true;
916  if ( mark != largeInteger ) {
917  //will be updated
918  work[mark] = value;
919  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
920  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
921 
922  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
923  added--;
924  } else {
925  thisPivotValue = value;
926  }
927  }
928  CoinBigIndex i;
929  for ( i = startColumn + 1; i < endColumn; i++ ) {
930  iRow = indexRowU[i];
931  value = elementU[i];
932  int mark = markRow[iRow];
933 
934  if ( mark == largeInteger+1 ) {
935  //keep
936  indexRowU[put] = iRow;
937  elementU[put] = value;
938  if ( checkLargest ) {
939  double absValue = fabs ( value );
940 
941  if ( absValue > largest ) {
942  largest = absValue;
943  positionLargest = put;
944  }
945  }
946  put++;
947  } else if ( mark != largeInteger ) {
948  //will be updated
949  work[mark] = value;
950  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
951  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
952 
953  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
954  added--;
955  } else {
956  thisPivotValue = value;
957  }
958  }
959  //slot in pivot
960  elementU[put] = elementU[startColumn];
961  indexRowU[put] = indexRowU[startColumn];
962  if ( positionLargest == startColumn ) {
963  positionLargest = put; //follow if was largest
964  }
965  put++;
966  elementU[startColumn] = thisPivotValue;
967  indexRowU[startColumn] = pivotRow;
968  //clean up counts
969  startColumn++;
970  numberInColumn[iColumn] = put - startColumn;
971  int * numberInColumnPlus = numberInColumnPlus_.array();
972  numberInColumnPlus[iColumn]++;
973  startColumnU[iColumn]++;
974  //how much space have we got
975  int next = nextColumn[iColumn];
976  CoinBigIndex space;
977 
978  space = startColumnU[next] - put - numberInColumnPlus[next];
979  //assume no zero elements
980  if ( numberInPivotColumn > space ) {
981  //getColumnSpace also moves fixed part
982  if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
983  return false;
984  }
985  //redo starts
986  if (positionLargest >= 0)
987  positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
988  startColumn = startColumnU[iColumn];
989  put = startColumn + numberInColumn[iColumn];
990  }
991  double tolerance = zeroTolerance_;
992 
993  int *nextCount = nextCount_.array();
994  for ( j = 0; j < numberInPivotColumn; j++ ) {
995  value = work[j] - thisPivotValue * multipliersL[j];
996  double absValue = fabs ( value );
997 
998  if ( absValue > tolerance ) {
999  work[j] = 0.0;
1000  assert (put<lengthAreaU_);
1001  elementU[put] = value;
1002  indexRowU[put] = indexL[j];
1003  if ( absValue > largest ) {
1004  largest = absValue;
1005  positionLargest = put;
1006  }
1007  put++;
1008  } else {
1009  work[j] = 0.0;
1010  added--;
1011  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1012  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1013 
1014  if ( temp2[word] & ( 1 << bit ) ) {
1015  //take out of row list
1016  iRow = indexL[j];
1017  CoinBigIndex start = startRowU[iRow];
1018  CoinBigIndex end = start + numberInRow[iRow];
1019  CoinBigIndex where = start;
1020 
1021  while ( indexColumnU[where] != iColumn ) {
1022  where++;
1023  } /* endwhile */
1024 #if DEBUG_COIN
1025  if ( where >= end ) {
1026  abort ( );
1027  }
1028 #endif
1029  indexColumnU[where] = indexColumnU[end - 1];
1030  numberInRow[iRow]--;
1031  } else {
1032  //make sure won't be added
1033  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1034  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1035 
1036  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1037  }
1038  }
1039  }
1040  numberInColumn[iColumn] = put - startColumn;
1041  //move largest
1042  if ( positionLargest >= 0 ) {
1043  value = elementU[positionLargest];
1044  iRow = indexRowU[positionLargest];
1045  elementU[positionLargest] = elementU[startColumn];
1046  indexRowU[positionLargest] = indexRowU[startColumn];
1047  elementU[startColumn] = value;
1048  indexRowU[startColumn] = iRow;
1049  }
1050  //linked list for column
1051  if ( nextCount[iColumn + numberRows_] != -2 ) {
1052  //modify linked list
1053  deleteLink ( iColumn + numberRows_ );
1054  addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1055  }
1056  temp2 += increment2;
1057  }
1058  //get space for row list
1059  unsigned int *putBase = workArea2;
1060  int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1061  int i = 0;
1062 
1063  // do linked lists and update counts
1064  while ( bigLoops ) {
1065  bigLoops--;
1066  int bit;
1067  for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
1068  unsigned int *putThis = putBase;
1069  int iRow = indexL[i];
1070 
1071  //get space
1072  int number = 0;
1073  int jColumn;
1074 
1075  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1076  unsigned int test = *putThis;
1077 
1078  putThis += increment2;
1079  test = 1 - ( ( test >> bit ) & 1 );
1080  number += test;
1081  }
1082  int next = nextRow[iRow];
1083  CoinBigIndex space;
1084 
1085  space = startRowU[next] - startRowU[iRow];
1086  number += numberInRow[iRow];
1087  if ( space < number ) {
1088  if ( !getRowSpace ( iRow, number ) ) {
1089  return false;
1090  }
1091  }
1092  // now do
1093  putThis = putBase;
1094  next = nextRow[iRow];
1095  number = numberInRow[iRow];
1096  CoinBigIndex end = startRowU[iRow] + number;
1097  int saveIndex = indexColumnU[startRowU[next]];
1098 
1099  //add in
1100  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1101  unsigned int test = *putThis;
1102 
1103  putThis += increment2;
1104  test = 1 - ( ( test >> bit ) & 1 );
1105  indexColumnU[end] = saveColumn[jColumn];
1106  end += test;
1107  }
1108  //put back next one in case zapped
1109  indexColumnU[startRowU[next]] = saveIndex;
1110  markRow[iRow] = static_cast<T>(largeInteger+1);
1111  number = end - startRowU[iRow];
1112  numberInRow[iRow] = number;
1113  deleteLink ( iRow );
1114  addLink ( iRow, number );
1115  }
1116  putBase++;
1117  } /* endwhile */
1118  int bit;
1119 
1120  for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
1121  unsigned int *putThis = putBase;
1122  int iRow = indexL[i];
1123 
1124  //get space
1125  int number = 0;
1126  int jColumn;
1127 
1128  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1129  unsigned int test = *putThis;
1130 
1131  putThis += increment2;
1132  test = 1 - ( ( test >> bit ) & 1 );
1133  number += test;
1134  }
1135  int next = nextRow[iRow];
1136  CoinBigIndex space;
1137 
1138  space = startRowU[next] - startRowU[iRow];
1139  number += numberInRow[iRow];
1140  if ( space < number ) {
1141  if ( !getRowSpace ( iRow, number ) ) {
1142  return false;
1143  }
1144  }
1145  // now do
1146  putThis = putBase;
1147  next = nextRow[iRow];
1148  number = numberInRow[iRow];
1149  CoinBigIndex end = startRowU[iRow] + number;
1150  int saveIndex;
1151 
1152  saveIndex = indexColumnU[startRowU[next]];
1153 
1154  //add in
1155  for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1156  unsigned int test = *putThis;
1157 
1158  putThis += increment2;
1159  test = 1 - ( ( test >> bit ) & 1 );
1160 
1161  indexColumnU[end] = saveColumn[jColumn];
1162  end += test;
1163  }
1164  indexColumnU[startRowU[next]] = saveIndex;
1165  markRow[iRow] = static_cast<T>(largeInteger+1);
1166  number = end - startRowU[iRow];
1167  numberInRow[iRow] = number;
1168  deleteLink ( iRow );
1169  addLink ( iRow, number );
1170  }
1171  markRow[pivotRow] = static_cast<T>(largeInteger+1);
1172  //modify linked list for pivots
1173  deleteLink ( pivotRow );
1174  deleteLink ( pivotColumn + numberRows_ );
1175  totalElements_ += added;
1176  return true;
1177 }
1178 
1179  /********************************* END LARGE TEMPLATE ********/
1181 protected:
1183 
1186  double pivotTolerance_;
1190 #ifndef COIN_FAST_CODE
1191  double slackValue_;
1193 #else
1194 #ifndef slackValue_
1195 #define slackValue_ -1.0
1196 #endif
1197 #endif
1198  double areaFactor_;
1201  double relaxCheck_;
1236  int status_;
1237 
1242  //int increasingRows_;
1243 
1248 
1251 
1254 
1257 
1261 
1264 
1267 
1270 
1273 
1276 
1279 
1282 
1285 
1288 
1291 
1294 
1297 
1300 
1303 
1306 
1309 
1311  //int baseU_;
1312 
1315 
1318 
1321 
1324 
1327 
1330 
1333 
1336 
1339 
1342 
1345 
1348 
1351 
1354 
1357 
1360 
1363 
1366 
1369 
1372 
1374  double * denseArea_;
1375 
1378 
1381 
1384 
1387 
1390 
1393 
1395  mutable double ftranCountInput_;
1396  mutable double ftranCountAfterL_;
1397  mutable double ftranCountAfterR_;
1398  mutable double ftranCountAfterU_;
1399  mutable double btranCountInput_;
1400  mutable double btranCountAfterU_;
1401  mutable double btranCountAfterR_;
1402  mutable double btranCountAfterL_;
1403 
1405  mutable int numberFtranCounts_;
1406  mutable int numberBtranCounts_;
1407 
1415 
1417  mutable bool collectStatistics_;
1418 
1421 
1424 
1427 
1430 
1433 
1439  int biasLU_;
1447 };
1448 // Dense coding
1449 #ifdef COIN_HAS_LAPACK
1450 #define DENSE_CODE 1
1451 /* Type of Fortran integer translated into C */
1452 #ifndef ipfint
1453 //typedef ipfint FORTRAN_INTEGER_TYPE ;
1454 typedef int ipfint;
1455 typedef const int cipfint;
1456 #endif
1457 #endif
1458 #endif
1459 // Extra for ugly include
1460 #ifdef UGLY_COIN_FACTOR_CODING
1461 #define FAC_UNSET (FAC_SET+1)
1462 {
1463  goodPivot=false;
1464  //store pivot columns (so can easily compress)
1465  CoinBigIndex startColumnThis = startColumn[iPivotColumn];
1466  CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
1467  int put = 0;
1468  CoinBigIndex startRowThis = startRow[iPivotRow];
1469  CoinBigIndex endRow = startRowThis + numberDoRow + 1;
1470  if ( pivotColumnPosition < 0 ) {
1471  for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1472  int iColumn = indexColumn[pivotColumnPosition];
1473  if ( iColumn != iPivotColumn ) {
1474  saveColumn[put++] = iColumn;
1475  } else {
1476  break;
1477  }
1478  }
1479  } else {
1480  for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
1481  saveColumn[put++] = indexColumn[i];
1482  }
1483  }
1484  assert (pivotColumnPosition<endRow);
1485  assert (indexColumn[pivotColumnPosition]==iPivotColumn);
1486  pivotColumnPosition++;
1487  for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1488  saveColumn[put++] = indexColumn[pivotColumnPosition];
1489  }
1490  //take out this bit of indexColumn
1491  int next = nextRow[iPivotRow];
1492  int last = lastRow[iPivotRow];
1493 
1494  nextRow[last] = next;
1495  lastRow[next] = last;
1496  nextRow[iPivotRow] = numberGoodU_; //use for permute
1497  lastRow[iPivotRow] = -2;
1498  numberInRow[iPivotRow] = 0;
1499  //store column in L, compress in U and take column out
1500  CoinBigIndex l = lengthL_;
1501  // **** HORRID coding coming up but a goto seems best!
1502  {
1503  if ( l + numberDoColumn > lengthAreaL_ ) {
1504  //need more memory
1505  if ((messageLevel_&4)!=0)
1506  printf("more memory needed in middle of invert\n");
1507  goto BAD_PIVOT;
1508  }
1509  //l+=currentAreaL_->elementByColumn-elementL;
1510  CoinBigIndex lSave = l;
1511 
1512  CoinBigIndex * startColumnL = startColumnL_.array();
1513  startColumnL[numberGoodL_] = l; //for luck and first time
1514  numberGoodL_++;
1515  startColumnL[numberGoodL_] = l + numberDoColumn;
1516  lengthL_ += numberDoColumn;
1517  if ( pivotRowPosition < 0 ) {
1518  for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1519  int iRow = indexRow[pivotRowPosition];
1520  if ( iRow != iPivotRow ) {
1521  indexRowL[l] = iRow;
1522  elementL[l] = element[pivotRowPosition];
1523  markRow[iRow] = l - lSave;
1524  l++;
1525  //take out of row list
1526  CoinBigIndex start = startRow[iRow];
1527  CoinBigIndex end = start + numberInRow[iRow];
1528  CoinBigIndex where = start;
1529 
1530  while ( indexColumn[where] != iPivotColumn ) {
1531  where++;
1532  } /* endwhile */
1533 #if DEBUG_COIN
1534  if ( where >= end ) {
1535  abort ( );
1536  }
1537 #endif
1538  indexColumn[where] = indexColumn[end - 1];
1539  numberInRow[iRow]--;
1540  } else {
1541  break;
1542  }
1543  }
1544  } else {
1545  CoinBigIndex i;
1546 
1547  for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
1548  int iRow = indexRow[i];
1549 
1550  markRow[iRow] = l - lSave;
1551  indexRowL[l] = iRow;
1552  elementL[l] = element[i];
1553  l++;
1554  //take out of row list
1555  CoinBigIndex start = startRow[iRow];
1556  CoinBigIndex end = start + numberInRow[iRow];
1557  CoinBigIndex where = start;
1558 
1559  while ( indexColumn[where] != iPivotColumn ) {
1560  where++;
1561  } /* endwhile */
1562 #if DEBUG_COIN
1563  if ( where >= end ) {
1564  abort ( );
1565  }
1566 #endif
1567  indexColumn[where] = indexColumn[end - 1];
1568  numberInRow[iRow]--;
1569  assert (numberInRow[iRow]>=0);
1570  }
1571  }
1572  assert (pivotRowPosition<endColumn);
1573  assert (indexRow[pivotRowPosition]==iPivotRow);
1574  CoinFactorizationDouble pivotElement = element[pivotRowPosition];
1575  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1576 
1577  pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1578  pivotRowPosition++;
1579  for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1580  int iRow = indexRow[pivotRowPosition];
1581 
1582  markRow[iRow] = l - lSave;
1583  indexRowL[l] = iRow;
1584  elementL[l] = element[pivotRowPosition];
1585  l++;
1586  //take out of row list
1587  CoinBigIndex start = startRow[iRow];
1588  CoinBigIndex end = start + numberInRow[iRow];
1589  CoinBigIndex where = start;
1590 
1591  while ( indexColumn[where] != iPivotColumn ) {
1592  where++;
1593  } /* endwhile */
1594 #if DEBUG_COIN
1595  if ( where >= end ) {
1596  abort ( );
1597  }
1598 #endif
1599  indexColumn[where] = indexColumn[end - 1];
1600  numberInRow[iRow]--;
1601  assert (numberInRow[iRow]>=0);
1602  }
1603  markRow[iPivotRow] = FAC_SET;
1604  //compress pivot column (move pivot to front including saved)
1605  numberInColumn[iPivotColumn] = 0;
1606  //use end of L for temporary space
1607  int *indexL = &indexRowL[lSave];
1608  CoinFactorizationDouble *multipliersL = &elementL[lSave];
1609 
1610  //adjust
1611  int j;
1612 
1613  for ( j = 0; j < numberDoColumn; j++ ) {
1614  multipliersL[j] *= pivotMultiplier;
1615  }
1616  //zero out fill
1617  CoinBigIndex iErase;
1618  for ( iErase = 0; iErase < increment2 * numberDoRow;
1619  iErase++ ) {
1620  workArea2[iErase] = 0;
1621  }
1622  CoinBigIndex added = numberDoRow * numberDoColumn;
1623  unsigned int *temp2 = workArea2;
1624  int * nextColumn = nextColumn_.array();
1625 
1626  //pack down and move to work
1627  int jColumn;
1628  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1629  int iColumn = saveColumn[jColumn];
1630  CoinBigIndex startColumnThis = startColumn[iColumn];
1631  CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
1632  int iRow = indexRow[startColumnThis];
1633  CoinFactorizationDouble value = element[startColumnThis];
1634  double largest;
1635  CoinBigIndex put = startColumnThis;
1636  CoinBigIndex positionLargest = -1;
1637  CoinFactorizationDouble thisPivotValue = 0.0;
1638 
1639  //compress column and find largest not updated
1640  bool checkLargest;
1641  int mark = markRow[iRow];
1642 
1643  if ( mark == FAC_UNSET ) {
1644  largest = fabs ( value );
1645  positionLargest = put;
1646  put++;
1647  checkLargest = false;
1648  } else {
1649  //need to find largest
1650  largest = 0.0;
1651  checkLargest = true;
1652  if ( mark != FAC_SET ) {
1653  //will be updated
1654  workArea[mark] = value;
1655  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1656  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1657 
1658  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1659  added--;
1660  } else {
1661  thisPivotValue = value;
1662  }
1663  }
1664  CoinBigIndex i;
1665  for ( i = startColumnThis + 1; i < endColumn; i++ ) {
1666  iRow = indexRow[i];
1667  value = element[i];
1668  int mark = markRow[iRow];
1669 
1670  if ( mark == FAC_UNSET ) {
1671  //keep
1672  indexRow[put] = iRow;
1673  element[put] = value;
1674  if ( checkLargest ) {
1675  double absValue = fabs ( value );
1676 
1677  if ( absValue > largest ) {
1678  largest = absValue;
1679  positionLargest = put;
1680  }
1681  }
1682  put++;
1683  } else if ( mark != FAC_SET ) {
1684  //will be updated
1685  workArea[mark] = value;
1686  int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1687  int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1688 
1689  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1690  added--;
1691  } else {
1692  thisPivotValue = value;
1693  }
1694  }
1695  //slot in pivot
1696  element[put] = element[startColumnThis];
1697  indexRow[put] = indexRow[startColumnThis];
1698  if ( positionLargest == startColumnThis ) {
1699  positionLargest = put; //follow if was largest
1700  }
1701  put++;
1702  element[startColumnThis] = thisPivotValue;
1703  indexRow[startColumnThis] = iPivotRow;
1704  //clean up counts
1705  startColumnThis++;
1706  numberInColumn[iColumn] = put - startColumnThis;
1707  int * numberInColumnPlus = numberInColumnPlus_.array();
1708  numberInColumnPlus[iColumn]++;
1709  startColumn[iColumn]++;
1710  //how much space have we got
1711  int next = nextColumn[iColumn];
1712  CoinBigIndex space;
1713 
1714  space = startColumn[next] - put - numberInColumnPlus[next];
1715  //assume no zero elements
1716  if ( numberDoColumn > space ) {
1717  //getColumnSpace also moves fixed part
1718  if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
1719  goto BAD_PIVOT;
1720  }
1721  //redo starts
1722  positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
1723  startColumnThis = startColumn[iColumn];
1724  put = startColumnThis + numberInColumn[iColumn];
1725  }
1726  double tolerance = zeroTolerance_;
1727 
1728  int *nextCount = nextCount_.array();
1729  for ( j = 0; j < numberDoColumn; j++ ) {
1730  value = workArea[j] - thisPivotValue * multipliersL[j];
1731  double absValue = fabs ( value );
1732 
1733  if ( absValue > tolerance ) {
1734  workArea[j] = 0.0;
1735  element[put] = value;
1736  indexRow[put] = indexL[j];
1737  if ( absValue > largest ) {
1738  largest = absValue;
1739  positionLargest = put;
1740  }
1741  put++;
1742  } else {
1743  workArea[j] = 0.0;
1744  added--;
1745  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1746  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1747 
1748  if ( temp2[word] & ( 1 << bit ) ) {
1749  //take out of row list
1750  iRow = indexL[j];
1751  CoinBigIndex start = startRow[iRow];
1752  CoinBigIndex end = start + numberInRow[iRow];
1753  CoinBigIndex where = start;
1754 
1755  while ( indexColumn[where] != iColumn ) {
1756  where++;
1757  } /* endwhile */
1758 #if DEBUG_COIN
1759  if ( where >= end ) {
1760  abort ( );
1761  }
1762 #endif
1763  indexColumn[where] = indexColumn[end - 1];
1764  numberInRow[iRow]--;
1765  } else {
1766  //make sure won't be added
1767  int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1768  int bit = j & COINFACTORIZATION_MASK_PER_INT;
1769 
1770  temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1771  }
1772  }
1773  }
1774  numberInColumn[iColumn] = put - startColumnThis;
1775  //move largest
1776  if ( positionLargest >= 0 ) {
1777  value = element[positionLargest];
1778  iRow = indexRow[positionLargest];
1779  element[positionLargest] = element[startColumnThis];
1780  indexRow[positionLargest] = indexRow[startColumnThis];
1781  element[startColumnThis] = value;
1782  indexRow[startColumnThis] = iRow;
1783  }
1784  //linked list for column
1785  if ( nextCount[iColumn + numberRows_] != -2 ) {
1786  //modify linked list
1787  deleteLink ( iColumn + numberRows_ );
1788  addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1789  }
1790  temp2 += increment2;
1791  }
1792  //get space for row list
1793  unsigned int *putBase = workArea2;
1794  int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1795  int i = 0;
1796 
1797  // do linked lists and update counts
1798  while ( bigLoops ) {
1799  bigLoops--;
1800  int bit;
1801  for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
1802  unsigned int *putThis = putBase;
1803  int iRow = indexL[i];
1804 
1805  //get space
1806  int number = 0;
1807  int jColumn;
1808 
1809  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1810  unsigned int test = *putThis;
1811 
1812  putThis += increment2;
1813  test = 1 - ( ( test >> bit ) & 1 );
1814  number += test;
1815  }
1816  int next = nextRow[iRow];
1817  CoinBigIndex space;
1818 
1819  space = startRow[next] - startRow[iRow];
1820  number += numberInRow[iRow];
1821  if ( space < number ) {
1822  if ( !getRowSpace ( iRow, number ) ) {
1823  goto BAD_PIVOT;
1824  }
1825  }
1826  // now do
1827  putThis = putBase;
1828  next = nextRow[iRow];
1829  number = numberInRow[iRow];
1830  CoinBigIndex end = startRow[iRow] + number;
1831  int saveIndex = indexColumn[startRow[next]];
1832 
1833  //add in
1834  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1835  unsigned int test = *putThis;
1836 
1837  putThis += increment2;
1838  test = 1 - ( ( test >> bit ) & 1 );
1839  indexColumn[end] = saveColumn[jColumn];
1840  end += test;
1841  }
1842  //put back next one in case zapped
1843  indexColumn[startRow[next]] = saveIndex;
1844  markRow[iRow] = FAC_UNSET;
1845  number = end - startRow[iRow];
1846  numberInRow[iRow] = number;
1847  deleteLink ( iRow );
1848  addLink ( iRow, number );
1849  }
1850  putBase++;
1851  } /* endwhile */
1852  int bit;
1853 
1854  for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
1855  unsigned int *putThis = putBase;
1856  int iRow = indexL[i];
1857 
1858  //get space
1859  int number = 0;
1860  int jColumn;
1861 
1862  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1863  unsigned int test = *putThis;
1864 
1865  putThis += increment2;
1866  test = 1 - ( ( test >> bit ) & 1 );
1867  number += test;
1868  }
1869  int next = nextRow[iRow];
1870  CoinBigIndex space;
1871 
1872  space = startRow[next] - startRow[iRow];
1873  number += numberInRow[iRow];
1874  if ( space < number ) {
1875  if ( !getRowSpace ( iRow, number ) ) {
1876  goto BAD_PIVOT;
1877  }
1878  }
1879  // now do
1880  putThis = putBase;
1881  next = nextRow[iRow];
1882  number = numberInRow[iRow];
1883  CoinBigIndex end = startRow[iRow] + number;
1884  int saveIndex;
1885 
1886  saveIndex = indexColumn[startRow[next]];
1887 
1888  //add in
1889  for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1890  unsigned int test = *putThis;
1891 
1892  putThis += increment2;
1893  test = 1 - ( ( test >> bit ) & 1 );
1894 
1895  indexColumn[end] = saveColumn[jColumn];
1896  end += test;
1897  }
1898  indexColumn[startRow[next]] = saveIndex;
1899  markRow[iRow] = FAC_UNSET;
1900  number = end - startRow[iRow];
1901  numberInRow[iRow] = number;
1902  deleteLink ( iRow );
1903  addLink ( iRow, number );
1904  }
1905  markRow[iPivotRow] = FAC_UNSET;
1906  //modify linked list for pivots
1907  deleteLink ( iPivotRow );
1908  deleteLink ( iPivotColumn + numberRows_ );
1909  totalElements_ += added;
1910  goodPivot= true;
1911  // **** UGLY UGLY UGLY
1912  }
1913  BAD_PIVOT:
1914 
1915  ;
1916 }
1917 #undef FAC_UNSET
1918 #endif
CoinBigIndex * array() const
Get Array.
CoinBigIndexArrayWithLength startRowU_
Start of each Row as pointer.
CoinBigIndex numberL_
Number in L.
int numberR_
Number in R.
CoinBigIndex factorElements_
Number of elements after factorization.
#define COINFACTORIZATION_BITS_PER_INT
CoinFactorizationDouble * elementU() const
Elements of U.
bool getColumnSpaceIterateR(int iColumn, double value, int iRow)
getColumnSpaceIterateR.
CoinIntArrayWithLength nextColumn_
Next Column in memory order.
bool pivotOneOtherRow(int pivotRow, int pivotColumn)
Pivots when just one other row so faster?
CoinBigIndex * startRowL() const
Start of each row in L.
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.
bool doForrestTomlin_
true if Forrest Tomlin update, false if PFI
CoinFactorizationDouble * pivotRegion() const
Returns address of pivot region.
int addColumn(CoinBigIndex numberElements, int indicesRow[], double elements[])
Adds one Column to basis, can increase size of basis.
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.
bool collectStatistics_
For statistics.
int status() const
Returns status.
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.
void deleteLink(int index)
Deletes a link in chain of equal counts.
void checkSparse()
See if worth going sparse.
int * pivotColumn() const
Returns address of pivotColumn region (also used for permuting)
CoinBigIndex lengthAreaL() const
Returns length of L area.
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.
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.
CoinBigIndex * startColumnL() const
Start of each column in L.
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 >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 >=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)
CoinFactorization()
Default constructor.
int deleteRow(int Row)
Deletes one Row from basis, returns rank.
CoinBigIndexArrayWithLength startColumnR_
Start of columns for R.
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.
double * denseArea_
Dense area.
CoinIntArrayWithLength markRow_
Marks rows to be updated.
int numberRowsExtra() const
Number of Rows after iterating.
void gutsOfInitialize(int type)
1 bit - tolerances etc, 2 more, 4 dummy arrays
bool pivot(int pivotRow, int pivotColumn, CoinBigIndex pivotRowPosition, CoinBigIndex pivotColumnPosition, CoinFactorizationDouble work[], unsigned int workArea2[], int increment2, T markRow[], int largeInteger)
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...
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.
double relaxCheck_
Relax check on accuracy in replaceColumn.
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.
CoinBigIndex lengthR_
Length of R stuff.
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 * 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 * 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.
CoinBigIndex * startColumnU() const
Start of each column in U.
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...
CoinBigIndex maximumU_
Maximum space used in U.
void updateColumnTransposeLSparse(CoinIndexedVector *region) const
Updates part of column transpose (BTRANL) when sparse (by Row)
CoinBigIndex lengthU_
Base of U is always 0.
void updateColumnPFI(CoinIndexedVector *regionSparse) const
Updates part of column PFI (FTRAN) (after rest)
double pivotTolerance_
Pivot tolerance.
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 <0 error, 0= finished.
int numberRowsExtra_
Number of Rows after iterating.
CoinBigIndex lengthL_
Length of L.
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 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...
double CoinFactorizationDouble
Definition: CoinTypes.hpp:54
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.
CoinBigIndex numberElementsL() const
Returns number in L area.
double getAccuracyCheck() const
void setBiasLU(int value)
int pivots() const
Returns number of pivots since factorization.
CoinBigIndex * version.
Indexed Vector.
CoinIntArrayWithLength nextCount_
Next Row/Column with count.
CoinBigIndexArrayWithLength startRowL_
Start of each row in L.
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.
CoinBigIndexArrayWithLength startColumnL_
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 denseThreshold_
Dense threshold.
void almostDestructor()
Delete all stuff (leaves as after CoinFactorization())
friend void CoinFactorizationUnitTest(const std::string &mpsDir)
#define COINFACTORIZATION_MASK_PER_INT
CoinIntArrayWithLength firstCount_
First Row/Column with count of k, can tell which by offset - Rows then Columns.
#define COIN_RESTRICT
void setNumberElementsU(CoinBigIndex value)
Setss number in U area.
CoinIntArrayWithLength sparse_
Sparse regions.
CoinFactorizationDouble * elementR_
Elements of R.
CoinBigIndex lengthAreaU_
Length of area reserved for U.
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.
int * indexRowL() const
Row indices of L.
void setStatus(int value)
Sets status.
CoinBigIndex getColumnSpaceIterate(int iColumn, double value, int iRow)
getColumnSpaceIterate.
CoinBigIndex lengthAreaR_
length of area reserved for R
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 addRow(CoinBigIndex numberElements, int indicesColumn[], double elements[])
Adds one Row to basis, can increase size of basis.
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)
CoinBigIndex lengthAreaU() const
Returns length of U area.
int maximumPivots() const
Maximum number of pivots between factorizations.
double adjustedAreaFactor() const
Returns areaFactor but adjusted for dense.
int factorSparseLarge()
Does sparse phase of factorization (for larger problems) return code is <0 error, 0= finished...
CoinBigIndex numberL() const
Number in L.
double ftranCountInput_
Below are all to collect.
int numberTrials_
0 - no increasing rows - no permutations, 1 - no increasing rows but permutations 2 - increasing rows...
CoinBigIndexArrayWithLength convertRowToColumnU_
Converts rows to columns in U.
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)
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.
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 <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.
CoinBigIndexArrayWithLength startColumnU_
Start of each column in U.
CoinBigIndex baseL_
Base of L.
CoinBigIndex numberElementsU() const
Returns number in U area.
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.
int CoinBigIndex
CoinIntArrayWithLength indexRowL_
Row indices of L.
CoinBigIndex numberCompressions_
Number of compressions done.
CoinBigIndex lengthAreaL_
Length of area reserved for L.
CoinIntArrayWithLength permute_
Permutation vector for pivot row order.
int add(CoinBigIndex numberElements, int indicesRow[], int indicesColumn[], double elements[])
Adds given elements to Basis and updates factorization, can increase size of basis.
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.
CoinBigIndex totalElements_
Number of elements in U (to go) or while iterating total overall.
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.
void getAreas(int numberRows, int numberColumns, CoinBigIndex maximumL, CoinBigIndex maximumU)
Gets space for a factorization, called by constructors.
CoinIntArrayWithLength numberInRow_
Number in each Row.
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.
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 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 "sparsish".
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.
void updateColumnLDensish(CoinIndexedVector *region, int *indexIn) const
Updates part of column (FTRANL) when densish.
CoinBigIndex numberCompressions() const
Number of compressions done.
CoinBigIndex baseL() const
Base of L.
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 <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)...
CoinBigIndex numberElementsR() const
Returns number in R area.
int * permuteBack() const
Returns address of permuteBack region.