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