1 /* $Id: CoinSimpFactorization.hpp 1448 2011-06-19 15:34:41Z stefan $ */
2 // Copyright (C) 2008, International Business Machines
3 // Corporation and others.  All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5 
6 /*
7    This is a simple factorization of the LP Basis
8  */
9 #ifndef CoinSimpFactorization_H
10 #define CoinSimpFactorization_H
11 
12 #include <iostream>
13 #include <string>
14 #include <cassert>
15 #include "CoinTypes.hpp"
16 #include "CoinIndexedVector.hpp"
17 #include "CoinDenseFactorization.hpp"
18 class CoinPackedMatrix;
19 
20 
21 /// pointers used during factorization
22 class FactorPointers{
23 public:
24     double *rowMax;
25     int *firstRowKnonzeros;
26     int *prevRow;
27     int *nextRow;
28     int *firstColKnonzeros;
29     int *prevColumn;
30     int *nextColumn;
31     int *newCols;
32     //constructor
33     FactorPointers( int numRows, int numCols, int *UrowLengths_, int *UcolLengths_ );
34     // destructor
35     ~ FactorPointers();
36 };
37 
38 class CoinSimpFactorization : public CoinOtherFactorization {
39    friend void CoinSimpFactorizationUnitTest( const std::string & mpsDir );
40 
41 public:
42 
43   /**@name Constructors and destructor and copy */
44   //@{
45   /// Default constructor
46   CoinSimpFactorization (  );
47   /// Copy constructor
48   CoinSimpFactorization ( const CoinSimpFactorization &other);
49 
50   /// Destructor
51   virtual ~CoinSimpFactorization (  );
52   /// = copy
53   CoinSimpFactorization & operator = ( const CoinSimpFactorization & other );
54   /// Clone
55   virtual CoinOtherFactorization * clone() const override ;
56   //@}
57 
58   /**@name Do factorization - public */
59   //@{
60   /// Gets space for a factorization
61   virtual void getAreas ( int numberRows,
62 		  int numberColumns,
63 		  CoinBigIndex maximumL,
64 		  CoinBigIndex maximumU ) override;
65 
66   /// PreProcesses column ordered copy of basis
67   virtual void preProcess ( ) override;
68   /** Does most of factorization returning status
69       0 - OK
70       -99 - needs more memory
71       -1 - singular - use numberGoodColumns and redo
72   */
73   virtual int factor ( ) override;
74   /// Does post processing on valid factorization - putting variables on correct rows
75   virtual void postProcess(const int * sequence, int * pivotVariable) override;
76   /// Makes a non-singular basis by replacing variables
77   virtual void makeNonSingular(int * sequence, int numberColumns) override;
78   //@}
79 
80   /**@name general stuff such as status */
81   //@{
82   /// Total number of elements in factorization
numberElements() const83   virtual inline int numberElements (  ) const override {
84     return numberRows_*(numberColumns_+numberPivots_);
85   }
86   /// Returns maximum absolute value in factorization
87   double maximumCoefficient() const;
88   //@}
89 
90   /**@name rank one updates which do exist */
91   //@{
92 
93   /** Replaces one Column to basis,
94    returns 0=OK, 1=Probably OK, 2=singular, 3=no room
95       If checkBeforeModifying is true will do all accuracy checks
96       before modifying factorization.  Whether to set this depends on
97       speed considerations.  You could just do this on first iteration
98       after factorization and thereafter re-factorize
99    partial update already in U */
100   virtual int replaceColumn ( CoinIndexedVector * regionSparse,
101 		      int pivotRow,
102 		      double pivotCheck ,
103 			      bool checkBeforeModifying=false,
104 			      double acceptablePivot=1.0e-8) override;
105   //@}
106 
107   /**@name various uses of factorization (return code number elements)
108    which user may want to know about */
109   //@{
110   /** Updates one column (FTRAN) from regionSparse2
111       Tries to do FT update
112       number returned is negative if no room
113       regionSparse starts as zero and is zero at end.
114       Note - if regionSparse2 packed on input - will be packed on output
115   */
116 
117     virtual int updateColumnFT ( CoinIndexedVector * regionSparse,
118 			 CoinIndexedVector * regionSparse2,
119 			 bool noPermute=false) override;
120 
121     /** This version has same effect as above with FTUpdate==false
122 	so number returned is always >=0 */
123     virtual int updateColumn ( CoinIndexedVector * regionSparse,
124 		       CoinIndexedVector * regionSparse2,
125 		       bool noPermute=false) const override;
126     /// does FTRAN on two columns
127     virtual int updateTwoColumnsFT(CoinIndexedVector * regionSparse1,
128 			   CoinIndexedVector * regionSparse2,
129 			   CoinIndexedVector * regionSparse3,
130 			   bool noPermute=false) override;
131     /// does updatecolumn if save==true keeps column for replace column
132     int upColumn ( CoinIndexedVector * regionSparse,
133 		   CoinIndexedVector * regionSparse2,
134 		   bool noPermute=false, bool save=false) const;
135     /** Updates one column (BTRAN) from regionSparse2
136 	regionSparse starts as zero and is zero at end
137 	Note - if regionSparse2 packed on input - will be packed on output
138     */
139     virtual int updateColumnTranspose ( CoinIndexedVector * regionSparse,
140 				CoinIndexedVector * regionSparse2) const override;
141     /// does updateColumnTranspose, the other is a wrapper
142     int upColumnTranspose ( CoinIndexedVector * regionSparse,
143 				CoinIndexedVector * regionSparse2) const;
144     //@}
145     /// *** Below this user may not want to know about
146 
147   /**@name various uses of factorization
148    which user may not want to know about (left over from my LP code) */
149   //@{
150   /// Get rid of all memory
clearArrays()151   inline void clearArrays() override
152   { gutsOfDestructor();}
153   /// Returns array to put basis indices in
indices() const154   inline int * indices() const override
155   { return reinterpret_cast<int *> (elements_+numberRows_*numberRows_);}
156   /// Returns permute in
permute() const157   virtual inline int * permute() const override
158   { return pivotRow_;}
159   //@}
160 
161   /// The real work of destructor
162   void gutsOfDestructor();
163   /// The real work of constructor
164   void gutsOfInitialize();
165   /// The real work of copy
166   void gutsOfCopy(const CoinSimpFactorization &other);
167 
168 
169     /// calls factorization
170   void factorize(int numberOfRows,
171 		 int numberOfColumns,
172 		 const int colStarts[],
173 		 const int indicesRow[],
174 		 const double elements[]);
175     /// main loop of factorization
176     int mainLoopFactor (FactorPointers &pointers );
177     /// copies L by rows
178     void copyLbyRows();
179     /// copies U by columns
180     void copyUbyColumns();
181     /// finds a pivot element using Markowitz count
182     int findPivot(FactorPointers &pointers, int &r, int &s, bool &ifSlack);
183     /// finds a pivot in a shortest column
184     int findPivotShCol(FactorPointers &pointers, int &r, int &s);
185     /// finds a pivot in the first column available
186     int findPivotSimp(FactorPointers &pointers, int &r, int &s);
187     /// does Gauss elimination
188     void GaussEliminate(FactorPointers &pointers, int &r, int &s);
189     /// finds short row that intersects a given column
190     int findShortRow(const int column, const int length, int &minRow,
191 		     int &minRowLength, FactorPointers &pointers);
192     /// finds short column that intersects a given row
193     int findShortColumn(const int row, const int length, int &minCol,
194 			int &minColLength, FactorPointers &pointers);
195     /// finds maximum absolute value in a row
196     double findMaxInRrow(const int row, FactorPointers &pointers);
197     /// does pivoting
198     void pivoting(const int pivotRow, const int pivotColumn,
199 		  const double invPivot, FactorPointers &pointers);
200     /// part of pivoting
201     void updateCurrentRow(const int pivotRow, const int row,
202 			  const double multiplier, FactorPointers &pointers,
203 			  int &newNonZeros);
204     /// allocates more space for L
205     void increaseLsize();
206     /// allocates more space for a row of U
207     void increaseRowSize(const int row, const int newSize);
208     /// allocates more space for a column of U
209     void increaseColSize(const int column, const int newSize, const bool b);
210     /// allocates more space for rows of U
211     void enlargeUrow(const int numNewElements);
212     /// allocates more space for columns of U
213     void enlargeUcol(const int numNewElements, const bool b);
214     /// finds a given row in a column
215     int findInRow(const int row, const int column);
216     /// finds a given column in a row
217     int findInColumn(const int column, const int row);
218     /// declares a row inactive
219     void removeRowFromActSet(const int row, FactorPointers &pointers);
220     /// declares a column inactive
221     void removeColumnFromActSet(const int column, FactorPointers &pointers);
222     /// allocates space for U
223     void allocateSpaceForU();
224     /// allocates several working arrays
225     void allocateSomeArrays();
226     /// initializes some numbers
227     void initialSomeNumbers();
228     /// solves L x = b
229     void Lxeqb(double *b) const;
230     /// same as above but with two rhs
231     void Lxeqb2(double *b1, double *b2) const;
232     /// solves U x = b
233     void Uxeqb(double *b, double *sol) const;
234     /// same as above but with two rhs
235     void Uxeqb2(double *b1, double *sol1, double *sol2, double *b2) const;
236     /// solves x L = b
237     void xLeqb(double *b) const;
238     /// solves x U = b
239     void xUeqb(double *b, double *sol) const;
240     /// updates factorization after a Simplex iteration
241     int LUupdate(int newBasicCol);
242     /// creates a new eta vector
243     void newEta(int row, int numNewElements);
244     /// makes a copy of row permutations
245     void copyRowPermutations();
246     /// solves H x = b, where H is a product of eta matrices
247     void Hxeqb(double *b) const;
248     /// same as above but with two rhs
249     void Hxeqb2(double *b1, double *b2) const;
250     /// solves x H = b
251     void xHeqb(double *b) const;
252     /// does FTRAN
253     void ftran(double *b, double *sol, bool save) const;
254     /// same as above but with two columns
255     void ftran2(double *b1, double *sol1, double *b2, double *sol2) const;
256     /// does BTRAN
257     void btran(double *b, double *sol) const;
258    ///---------------------------------------
259 
260 
261 
262   //@}
263 protected:
264   /** Returns accuracy status of replaceColumn
265       returns 0=OK, 1=Probably OK, 2=singular */
266   int checkPivot(double saveFromU, double oldPivot) const;
267 ////////////////// data //////////////////
268 protected:
269 
270   /**@name data */
271   //@{
272     /// work array (should be initialized to zero)
273     double *denseVector_;
274     /// work array
275     double *workArea2_;
276     /// work array
277     double *workArea3_;
278     /// array of labels (should be initialized to zero)
279     int *vecLabels_;
280     /// array of indices
281     int *indVector_;
282 
283     /// auxiliary vector
284     double *auxVector_;
285     /// auxiliary vector
286     int *auxInd_;
287 
288     /// vector to keep for LUupdate
289     double *vecKeep_;
290     /// indices of this vector
291     int *indKeep_;
292     /// number of nonzeros
293     mutable int keepSize_;
294 
295 
296 
297     /// Starts of the rows of L
298     int *LrowStarts_;
299     /// Lengths of the rows of L
300     int *LrowLengths_;
301     /// L by rows
302     double *Lrows_;
303     /// indices in the rows of L
304     int *LrowInd_;
305     /// Size of Lrows_;
306     int LrowSize_;
307     /// Capacity of Lrows_
308     int LrowCap_;
309 
310     /// Starts of the columns of L
311     int *LcolStarts_;
312     /// Lengths of the columns of L
313     int *LcolLengths_;
314     /// L by columns
315     double *Lcolumns_;
316     /// indices in the columns of L
317     int *LcolInd_;
318     /// numbers of elements in L
319     int LcolSize_;
320     /// maximum capacity of L
321     int LcolCap_;
322 
323 
324     /// Starts of the rows of U
325     int *UrowStarts_;
326     /// Lengths of the rows of U
327     int *UrowLengths_;
328 #ifdef COIN_SIMP_CAPACITY
329     /// Capacities of the rows of U
330     int *UrowCapacities_;
331 #endif
332     /// U by rows
333     double *Urows_;
334     /// Indices in the rows of U
335     int *UrowInd_;
336     /// maximum capacity of Urows
337     int UrowMaxCap_;
338     /// number of used places in Urows
339     int UrowEnd_;
340     /// first row in U
341     int firstRowInU_;
342     /// last row in U
343     int lastRowInU_;
344     /// previous row in U
345     int *prevRowInU_;
346     /// next row in U
347     int *nextRowInU_;
348 
349     /// Starts of the columns of U
350     int *UcolStarts_;
351     /// Lengths of the columns of U
352     int *UcolLengths_;
353 #ifdef COIN_SIMP_CAPACITY
354     /// Capacities of the columns of U
355     int *UcolCapacities_;
356 #endif
357     /// U by columns
358     double *Ucolumns_;
359     /// Indices in the columns of U
360     int *UcolInd_;
361     /// previous column in U
362     int *prevColInU_;
363     /// next column in U
364     int *nextColInU_;
365     /// first column in U
366     int firstColInU_;
367     /// last column in U
368     int lastColInU_;
369     /// maximum capacity of Ucolumns_
370     int UcolMaxCap_;
371     /// last used position in Ucolumns_
372     int UcolEnd_;
373     /// indicator of slack variables
374     int *colSlack_;
375 
376     /// inverse values of the elements of diagonal of U
377     double *invOfPivots_;
378 
379     /// permutation of columns
380     int *colOfU_;
381     /// position of column after permutation
382     int *colPosition_;
383     /// permutations of rows
384     int *rowOfU_;
385     /// position of row after permutation
386     int *rowPosition_;
387     /// permutations of rows during LUupdate
388     int *secRowOfU_;
389     /// position of row after permutation during LUupdate
390     int *secRowPosition_;
391 
392     /// position of Eta vector
393     int *EtaPosition_;
394     /// Starts of eta vectors
395     int *EtaStarts_;
396     /// Lengths of eta vectors
397     int *EtaLengths_;
398     /// columns of eta vectors
399     int *EtaInd_;
400     /// elements of eta vectors
401     double *Eta_;
402     /// number of elements in Eta_
403     int EtaSize_;
404     /// last eta row
405     int lastEtaRow_;
406     /// maximum number of eta vectors
407     int maxEtaRows_;
408     /// Capacity of Eta_
409     int EtaMaxCap_;
410 
411     /// minimum storage increase
412     int minIncrease_;
413     /// maximum size for the diagonal of U after update
414     double updateTol_;
415     /// do Shul heuristic
416     bool doSuhlHeuristic_;
417     /// maximum of U
418     double maxU_;
419     /// bound on the growth rate
420     double maxGrowth_;
421     /// maximum of A
422     double maxA_;
423     /// maximum number of candidates for pivot
424     int pivotCandLimit_;
425     /// number of slacks in basis
426     int numberSlacks_;
427     /// number of slacks in irst basis
428     int firstNumberSlacks_;
429     //@}
430 };
431 #endif
432