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"
18class CoinPackedMatrix;
19
20
21/// pointers used during factorization
22class FactorPointers{
23public:
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
38class CoinSimpFactorization : public CoinOtherFactorization {
39 friend void CoinSimpFactorizationUnitTest( const std::string & mpsDir );
40
41public:
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
83 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
151 inline void clearArrays() override
152 { gutsOfDestructor();}
153 /// Returns array to put basis indices in
154 inline int * indices() const override
155 { return reinterpret_cast<int *> (elements_+numberRows_*numberRows_);}
156 /// Returns permute in
157 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 //@}
263protected:
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 //////////////////
268protected:
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