1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the class library */ 4 /* SoPlex --- the Sequential object-oriented simPlex. */ 5 /* */ 6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file clufactor.h 26 * @brief Implementation of sparse LU factorization. 27 */ 28 #ifndef _CLUFACTOR_H_ 29 #define _CLUFACTOR_H_ 30 31 #include "soplex/spxdefines.h" 32 #include "soplex/slinsolver.h" 33 #include "soplex/timer.h" 34 #include "soplex/svector.h" 35 36 #include "vector" 37 38 #define SOPLEX_WITH_L_ROWS 1 39 40 namespace soplex 41 { 42 /**@brief Implementation of sparse LU factorization. 43 * @ingroup Algo 44 * 45 * This class implements a sparse LU factorization with either 46 * FOREST-TOMLIN or ETA updates, using dynamic Markowitz pivoting. 47 */ 48 template <class R> 49 class CLUFactor 50 { 51 public: 52 53 /// set tolerances 54 void setTolerances(std::shared_ptr<Tolerances> tolerances) 55 { 56 this->_tolerances = tolerances; 57 } 58 /// get tolerances 59 const std::shared_ptr<Tolerances> tolerances() const 60 { 61 return _tolerances; 62 } 63 64 //---------------------------------------- 65 /**@name Public types */ 66 ///@{ 67 /** Doubly linked ring structure for garbage collection of column or 68 * row file in working matrix. 69 */ 70 struct Dring 71 { 72 Dring* next; 73 Dring* prev; 74 int idx; 75 }; 76 77 /// Pivot Ring 78 class Pring 79 { 80 public: 81 Pring* next; ///< 82 Pring* prev; ///< 83 int idx; ///< index of pivot row 84 int pos; ///< position of pivot column in row 85 int mkwtz; ///< markowitz number of pivot 86 87 Pring() : next(0), prev(0) ///< constructor 88 { 89 mkwtz = -1; 90 idx = -1; 91 pos = -1; 92 } 93 94 private: 95 Pring(const Pring&); ///< blocked copy constructor 96 Pring& operator= (const Pring&); ///< blocked assignment operator 97 }; 98 ///@} 99 100 protected: 101 102 //---------------------------------------- 103 /**@name Protected types */ 104 ///@{ 105 /// Temporary data structures. 106 class Temp 107 { 108 public: 109 int* s_mark; ///< marker 110 std::vector<R> s_max; ///< maximum absolute value per row (or -1) 111 int* s_cact; ///< lengths of columns of active submatrix 112 int stage; ///< stage of the structure 113 Pring pivots; ///< ring of selected pivot rows 114 Pring* pivot_col; ///< column index handlers for R linked list 115 Pring* pivot_colNZ; ///< lists for columns to number of nonzeros 116 Pring* pivot_row; ///< row index handlers for R linked list 117 Pring* pivot_rowNZ; ///< lists for rows to number of nonzeros 118 119 Temp(); ///< constructor 120 ~Temp(); ///< destructor 121 void init(int p_dim); ///< initialization 122 void clear(); ///< clears the structure 123 124 private: 125 Temp(const Temp&); ///< blocked copy constructor 126 Temp& operator= (const Temp&); ///< blocked assignment operator 127 }; 128 129 /// Data structures for saving the row and column permutations. 130 struct Perm 131 { 132 int* orig; ///< orig[p] original index from p 133 int* perm; ///< perm[i] permuted index from i 134 }; 135 136 /// Data structures for saving the working matrix and U factor. 137 struct U 138 { 139 /// 140 struct Row 141 { 142 Dring list; /*!< \brief Double linked ringlist of VectorBase<R> 143 indices in the order they appear 144 in the row file */ 145 Dring* elem; ///< %Array of ring elements. 146 int size; ///< size of arrays val and idx 147 int used; ///< used entries of arrays idx and val 148 std::vector<R> val; ///< hold nonzero values 149 int* idx; ///< hold column indices of nonzeros 150 int* start; ///< starting positions in val and idx 151 int* len; ///< used nonzeros per row vectors 152 int* max; /*!< \brief maximum available nonzeros per row: 153 start[i] + max[i] == start[elem[i].next->idx] 154 len[i] <= max[i]. */ 155 } row; 156 157 /// 158 struct Col 159 { 160 Dring list; /*!< \brief Double linked ringlist of VectorBase<R> 161 indices in the order they appear 162 in the column file */ 163 Dring* elem; ///< %Array of ring elements. 164 int size; ///< size of array idx 165 int used; ///< used entries of array idx 166 int* idx; ///< hold row indices of nonzeros 167 std::vector<R> val; /*!< \brief hold nonzero values: this is only initialized 168 in the end of the factorization with DEFAULT 169 updates. */ 170 int* start; ///< starting positions in val and idx 171 int* len; ///< used nonzeros per column vector 172 int* max; /*!< \brief maximum available nonzeros per colunn: 173 start[i] + max[i] == start[elem[i].next->idx] 174 len[i] <= max[i]. */ 175 } col; 176 }; 177 178 179 /// Data structures for saving the working matrix and L factor. 180 struct L 181 { 182 int size; ///< size of arrays val and idx 183 std::vector<R> val; ///< values of L vectors 184 int* idx; ///< indices of L vectors 185 int startSize; ///< size of array start 186 int firstUpdate; ///< number of first update L vector 187 int firstUnused; ///< number of first unused L vector 188 int* start; ///< starting positions in val and idx 189 int* row; ///< column indices of L vectors 190 int updateType; ///< type of updates to be used. 191 192 /* The following arrays have length |firstUpdate|, since they keep 193 * rows of the L-vectors occuring during the factorization (without 194 * updates), only: 195 */ 196 std::vector<R> rval; ///< values of rows of L 197 int* ridx; ///< indices of rows of L 198 int* rbeg; ///< start of rows in rval and ridx 199 int* rorig; ///< original row permutation 200 int* rperm; ///< original row permutation 201 }; 202 ///@} 203 204 //---------------------------------------- 205 /**@name Protected data */ 206 ///@{ 207 typename SLinSolver<R>::Status stat; ///< Status indicator. 208 209 int thedim; ///< dimension of factorized matrix 210 int nzCnt; ///< number of nonzeros in U 211 R initMaxabs; ///< maximum abs number in initail Matrix 212 R maxabs; ///< maximum abs number in L and U 213 214 R rowMemMult; ///< factor of minimum Memory * number of nonzeros 215 R colMemMult; ///< factor of minimum Memory * number of nonzeros 216 R lMemMult; ///< factor of minimum Memory * number of nonzeros 217 218 Perm row; ///< row permutation matrices 219 Perm col; ///< column permutation matrices 220 221 L l; ///< L matrix 222 std::vector<R> diag; ///< Array of pivot elements 223 U u; ///< U matrix 224 225 R* work; ///< Working array: must always be left as 0! 226 227 Timer* factorTime; ///< Time spent in factorizations 228 int factorCount; ///< Number of factorizations 229 int hugeValues; ///< number of times huge values occurred during solve (only used in debug mode) 230 std::shared_ptr<Tolerances> _tolerances; ///< Tolerances for the factorization 231 ///@} 232 233 private: 234 235 //---------------------------------------- 236 /**@name Private data */ 237 ///@{ 238 Temp temp; ///< Temporary storage 239 ///@} 240 241 //---------------------------------------- 242 /**@name Solving 243 These helper methods are used during the factorization process. 244 The solve*-methods solve lower and upper triangular systems from 245 the left or from the right, respectively The methods with '2' in 246 the end solve two systems at the same time. The methods with 247 "Eps" in the end consider elements smaller then the passed epsilon 248 as zero. 249 */ 250 ///@{ 251 /// 252 void solveUright(R* wrk, R* vec) const; 253 /// 254 int solveUrightEps(R* vec, int* nonz, R eps, R* rhs); 255 /// 256 void solveUright2(R* work1, R* vec1, R* work2, R* vec2); 257 /// 258 int solveUright2eps(R* work1, R* vec1, R* work2, R* vec2, int* nonz, R eps); 259 /// 260 void solveLright2(R* vec1, R* vec2); 261 /// 262 void solveUpdateRight(R* vec); 263 /// 264 void solveUpdateRight2(R* vec1, R* vec2); 265 /// 266 void solveUleft(R* work, R* vec); 267 /// 268 void solveUleft2(R* work1, R* vec1, R* work2, R* vec2); 269 /// 270 int solveLleft2forest(R* vec1, int* /* nonz */, R* vec2, R /* eps */); 271 /// 272 void solveLleft2(R* vec1, int* /* nonz */, R* vec2, R /* eps */); 273 /// 274 int solveLleftForest(R* vec, int* /* nonz */, R /* eps */); 275 /// 276 void solveLleft(R* vec) const; 277 /// 278 int solveLleftEps(R* vec, int* nonz, R eps); 279 /// 280 void solveUpdateLeft(R* vec); 281 /// 282 void solveUpdateLeft2(R* vec1, R* vec2); 283 284 void inline updateSolutionVectorLright(R change, int j, R& vec, int* idx, int& nnz); 285 /// 286 void vSolveLright(R* vec, int* ridx, int& rn, R eps); 287 /// 288 void vSolveLright2(R* vec, int* ridx, int& rn, R eps, 289 R* vec2, int* ridx2, int& rn2, R eps2); 290 /// 291 void vSolveLright3(R* vec, int* ridx, int& rn, R eps, 292 R* vec2, int* ridx2, int& rn2, R eps2, 293 R* vec3, int* ridx3, int& rn3, R eps3); 294 /// 295 int vSolveUright(R* vec, int* vidx, R* rhs, int* ridx, int rn, R eps); 296 /// 297 void vSolveUrightNoNZ(R* vec, R* rhs, int* ridx, int rn, R eps); 298 /// 299 int vSolveUright2(R* vec, int* vidx, R* rhs, int* ridx, int rn, R eps, 300 R* vec2, R* rhs2, int* ridx2, int rn2, R eps2); 301 /// 302 int vSolveUpdateRight(R* vec, int* ridx, int n, R eps); 303 /// 304 void vSolveUpdateRightNoNZ(R* vec, R /*eps*/); 305 /// 306 int solveUleft(R eps, R* vec, int* vecidx, R* rhs, int* rhsidx, int rhsn); 307 /// 308 void solveUleftNoNZ(R eps, R* vec, R* rhs, int* rhsidx, int rhsn); 309 /// 310 int solveLleftForest(R eps, R* vec, int* nonz, int n); 311 /// 312 void solveLleftForestNoNZ(R* vec); 313 /// 314 int solveLleft(R eps, R* vec, int* nonz, int rn); 315 /// 316 void solveLleftNoNZ(R* vec); 317 /// 318 int solveUpdateLeft(R eps, R* vec, int* nonz, int n); 319 320 /// 321 void forestPackColumns(); 322 /// 323 void forestMinColMem(int size); 324 /// 325 void forestReMaxCol(int col, int len); 326 327 /// 328 void initPerm(); 329 /// 330 void initFactorMatrix(const SVectorBase<R>** vec, const R eps); 331 /// 332 void minLMem(int size); 333 /// 334 void setPivot(const int p_stage, const int p_col, const int p_row, const R val); 335 /// 336 void colSingletons(); 337 /// 338 void rowSingletons(); 339 340 /// 341 void initFactorRings(); 342 /// 343 void freeFactorRings(); 344 345 /// 346 int setupColVals(); 347 /// 348 void setupRowVals(); 349 350 /// 351 void eliminateRowSingletons(); 352 /// 353 void eliminateColSingletons(); 354 /// 355 void selectPivots(R threshold); 356 /// 357 int updateRow(int r, int lv, int prow, int pcol, R pval, R eps); 358 359 /// 360 void eliminatePivot(int prow, int pos, R eps); 361 /// 362 void eliminateNucleus(const R eps, const R threshold); 363 /// 364 void minRowMem(int size); 365 /// 366 void minColMem(int size); 367 /// 368 void remaxCol(int p_col, int len); 369 /// 370 void packRows(); 371 /// 372 void packColumns(); 373 /// 374 void remaxRow(int p_row, int len); 375 /// 376 int makeLvec(int p_len, int p_row); 377 ///@} 378 379 protected: 380 381 //---------------------------------------- 382 /**@name Solver methods */ 383 ///@{ 384 /// 385 void solveLright(R* vec); 386 /// 387 int solveRight4update(R* vec, int* nonz, R eps, R* rhs, 388 R* forest, int* forestNum, int* forestIdx); 389 /// 390 void solveRight(R* vec, R* rhs); 391 /// 392 int solveRight2update(R* vec1, R* vec2, R* rhs1, 393 R* rhs2, int* nonz, R eps, R* forest, int* forestNum, int* forestIdx); 394 /// 395 void solveRight2(R* vec1, R* vec2, R* rhs1, R* rhs2); 396 /// 397 void solveLeft(R* vec, R* rhs); 398 /// 399 int solveLeftEps(R* vec, R* rhs, int* nonz, R eps); 400 /// 401 int solveLeft2(R* vec1, int* nonz, R* vec2, R eps, R* rhs1, R* rhs2); 402 403 /// 404 int vSolveRight4update(R eps, 405 R* vec, int* idx, /* result */ 406 R* rhs, int* ridx, int rn, /* rhs & Forest */ 407 R* forest, int* forestNum, int* forestIdx); 408 /// 409 int vSolveRight4update2(R eps, 410 R* vec, int* idx, /* result1 */ 411 R* rhs, int* ridx, int rn, /* rhs1 */ 412 R* vec2, R eps2, /* result2 */ 413 R* rhs2, int* ridx2, int rn2, /* rhs2 */ 414 R* forest, int* forestNum, int* forestIdx); 415 /// sparse version of above method 416 void vSolveRight4update2sparse( 417 R eps, R* vec, int* idx, /* result1 */ 418 R* rhs, int* ridx, int& rn, /* rhs1 */ 419 R eps2, R* vec2, int* idx2, /* result2 */ 420 R* rhs2, int* ridx2, int& rn2, /* rhs2 */ 421 R* forest, int* forestNum, int* forestIdx); 422 /// 423 int vSolveRight4update3(R eps, 424 R* vec, int* idx, /* result1 */ 425 R* rhs, int* ridx, int rn, /* rhs1 */ 426 R* vec2, R eps2, /* result2 */ 427 R* rhs2, int* ridx2, int rn2, /* rhs2 */ 428 R* vec3, R eps3, /* result3 */ 429 R* rhs3, int* ridx3, int rn3, /* rhs3 */ 430 R* forest, int* forestNum, int* forestIdx); 431 /// sparse version of above method 432 void vSolveRight4update3sparse( 433 R eps, R* vec, int* idx, /* result1 */ 434 R* rhs, int* ridx, int& rn, /* rhs1 */ 435 R eps2, R* vec2, int* idx2, /* result2 */ 436 R* rhs2, int* ridx2, int& rn2, /* rhs2 */ 437 R eps3, R* vec3, int* idx3, /* result3 */ 438 R* rhs3, int* ridx3, int& rn3, /* rhs3 */ 439 R* forest, int* forestNum, int* forestIdx); 440 /// 441 void vSolveRightNoNZ(R* vec, R eps, /* result */ 442 R* rhs, int* ridx, int rn); /* rhs */ 443 /// 444 int vSolveLeft(R eps, 445 R* vec, int* idx, /* result */ 446 R* rhs, int* ridx, int rn); /* rhs */ 447 /// 448 void vSolveLeftNoNZ(R eps, 449 R* vec, /* result */ 450 R* rhs, int* ridx, int rn); /* rhs */ 451 /// 452 int vSolveLeft2(R eps, 453 R* vec, int* idx, /* result */ 454 R* rhs, int* ridx, int rn, /* rhs */ 455 R* vec2, /* result2 */ 456 R* rhs2, int* ridx2, int rn2); /* rhs2 */ 457 /// sparse version of solving 2 systems of equations 458 void vSolveLeft2sparse(R eps, 459 R* vec, int* idx, /* result */ 460 R* rhs, int* ridx, int& rn, /* rhs */ 461 R* vec2, int* idx2, /* result2 */ 462 R* rhs2, int* ridx2, int& rn2); /* rhs2 */ 463 /// 464 int vSolveLeft3(R eps, 465 R* vec, int* idx, /* result */ 466 R* rhs, int* ridx, int rn, /* rhs */ 467 R* vec2, /* result2 */ 468 R* rhs2, int* ridx2, int rn2, /* rhs2 */ 469 R* vec3, /* result3 */ 470 R* rhs3, int* ridx3, int rn3); /* rhs3 */ 471 /// sparse version of solving 3 systems of equations 472 void vSolveLeft3sparse(R eps, 473 R* vec, int* idx, /* result */ 474 R* rhs, int* ridx, int& rn, /* rhs */ 475 R* vec2, int* idx2, /* result2 */ 476 R* rhs2, int* ridx2, int& rn2, /* rhs2 */ 477 R* vec3, int* idx3, /* result2 */ 478 R* rhs3, int* ridx3, int& rn3); /* rhs2 */ 479 480 void forestUpdate(int col, R* work, int num, int* nonz); 481 482 void update(int p_col, R* p_work, const int* p_idx, int num); 483 void updateNoClear(int p_col, const R* p_work, const int* p_idx, int num); 484 485 /// 486 void factor(const SVectorBase<R>** vec, ///< Array of column VectorBase<R> pointers 487 R threshold, ///< pivoting threshold 488 R eps); ///< epsilon for zero detection 489 ///@} 490 491 //---------------------------------------- 492 /**@name Debugging */ 493 ///@{ 494 /// 495 void dump() const; 496 497 /// 498 bool isConsistent() const; 499 ///@} 500 }; 501 502 } // namespace soplex 503 504 // For general templated functions 505 #include "clufactor.hpp" 506 507 #endif // _CLUFACTOR_H_ 508