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_rational.h 26 * @brief Implementation of sparse LU factorization with Rational precision. 27 */ 28 #ifndef _CLUFACTOR_RATIONAL_H_ 29 #define _CLUFACTOR_RATIONAL_H_ 30 31 #include "soplex/spxdefines.h" 32 #include "soplex/slinsolver_rational.h" 33 #include "soplex/timer.h" 34 #include "soplex/svector.h" 35 #include "soplex/rational.h" 36 #include "soplex/basevectors.h" 37 38 #define SOPLEX_WITH_L_ROWS 1 39 40 namespace soplex 41 { 42 /**@brief Implementation of sparse LU factorization with Rational precision. 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 class CLUFactorRational 49 { 50 public: 51 52 //---------------------------------------- 53 /**@name Public types */ 54 ///@{ 55 /** Doubly linked ring structure for garbage collection of column or 56 * row file in working matrix. 57 */ 58 struct Dring 59 { 60 Dring* next; 61 Dring* prev; 62 int idx; 63 }; 64 65 /// Pivot Ring 66 class Pring 67 { 68 public: 69 Pring* next; ///< 70 Pring* prev; ///< 71 int idx; ///< index of pivot row 72 int pos; ///< position of pivot column in row 73 int mkwtz; ///< markowitz number of pivot 74 75 Pring() ///< constructor 76 : next(0) 77 , prev(0) 78 , idx(0) 79 , pos(0) 80 , mkwtz(0) 81 {} 82 83 private: 84 Pring(const Pring&); ///< blocked copy constructor 85 Pring& operator= (const Pring&); ///< blocked assignment operator 86 }; 87 ///@} 88 89 protected: 90 91 //---------------------------------------- 92 /**@name Protected types */ 93 ///@{ 94 /// Temporary data structures. 95 class Temp 96 { 97 public: 98 int* s_mark; ///< marker 99 VectorRational s_max; ///< maximum absolute value per row (or -1) 100 int* s_cact; ///< lengths of columns of active submatrix 101 int stage; ///< stage of the structure 102 Pring pivots; ///< ring of selected pivot rows 103 Pring* pivot_col; ///< column index handlers for Real linked list 104 Pring* pivot_colNZ; ///< lists for columns to number of nonzeros 105 Pring* pivot_row; ///< row index handlers for Real linked list 106 Pring* pivot_rowNZ; ///< lists for rows to number of nonzeros 107 108 Temp(); ///< constructor 109 ~Temp(); ///< destructor 110 void init(int p_dim); ///< initialization 111 void clear(); ///< clears the structure 112 113 private: 114 Temp(const Temp&); ///< blocked copy constructor 115 Temp& operator= (const Temp&); ///< blocked assignment operator 116 }; 117 118 /// Data structures for saving the row and column permutations. 119 struct Perm 120 { 121 int* orig; ///< orig[p] original index from p 122 int* perm; ///< perm[i] permuted index from i 123 }; 124 125 /// Data structures for saving the working matrix and U factor. 126 struct U 127 { 128 /// 129 struct Row 130 { 131 Dring list; /*!< \brief Double linked ringlist of vector 132 indices in the order they appear 133 in the row file */ 134 Dring* elem; ///< %Array of ring elements. 135 int used; ///< used entries of arrays idx and val 136 VectorRational val; ///< hold nonzero values 137 int* idx; ///< array of length val.dim() to hold column indices of nonzeros in val 138 int* start; ///< starting positions in val and idx 139 int* len; ///< used nonzeros per row vectors 140 int* max; /*!< \brief maximum available nonzeros per row: 141 start[i] + max[i] == start[elem[i].next->idx] 142 len[i] <= max[i]. */ 143 } row; 144 145 /// 146 struct Col 147 { 148 Dring list; /*!< \brief Double linked ringlist of vector 149 indices in the order they appear 150 in the column file */ 151 Dring* elem; ///< %Array of ring elements. 152 int size; ///< size of array idx 153 int used; ///< used entries of array idx 154 int* idx; ///< hold row indices of nonzeros 155 VectorRational val; /*!< \brief hold nonzero values: this is only initialized 156 in the end of the factorization with DEFAULT 157 updates. */ 158 int* start; ///< starting positions in val and idx 159 int* len; ///< used nonzeros per column vector 160 int* max; /*!< \brief maximum available nonzeros per colunn: 161 start[i] + max[i] == start[elem[i].next->idx] 162 len[i] <= max[i]. */ 163 } col; 164 }; 165 166 167 /// Data structures for saving the working matrix and L factor. 168 struct L 169 { 170 VectorRational val; ///< values of L vectors 171 int* idx; ///< array of size val.dim() storing indices of L vectors 172 int startSize; ///< size of array start 173 int firstUpdate; ///< number of first update L vector 174 int firstUnused; ///< number of first unused L vector 175 int* start; ///< starting positions in val and idx 176 int* row; ///< column indices of L vectors 177 int updateType; ///< type of updates to be used. 178 179 /* The following arrays have length |firstUpdate|, since they keep 180 * rows of the L-vectors occuring during the factorization (without 181 * updates), only: 182 */ 183 VectorRational rval; ///< values of rows of L 184 int* ridx; ///< indices of rows of L 185 int* rbeg; ///< start of rows in rval and ridx 186 int* rorig; ///< original row permutation 187 int* rperm; ///< original row permutation 188 }; 189 ///@} 190 191 //---------------------------------------- 192 /**@name Protected data */ 193 ///@{ 194 SLinSolverRational::Status stat; ///< Status indicator. 195 196 int thedim; ///< dimension of factorized matrix 197 int nzCnt; ///< number of nonzeros in U 198 Rational initMaxabs; ///< maximum abs number in initail Matrix 199 Rational maxabs; ///< maximum abs number in L and U 200 201 Real rowMemMult; ///< factor of minimum Memory * number of nonzeros 202 Real colMemMult; ///< factor of minimum Memory * number of nonzeros 203 Real lMemMult; ///< factor of minimum Memory * number of nonzeros 204 205 Perm row; ///< row permutation matrices 206 Perm col; ///< column permutation matrices 207 208 L l; ///< L matrix 209 VectorRational diag; ///< Array of pivot elements 210 U u; ///< U matrix 211 212 Rational* work; ///< Working array: must always be left as 0! 213 214 Timer* factorTime; ///< Time spent in factorizations 215 int factorCount; ///< Number of factorizations 216 Real timeLimit; ///< Time limit on factorization or solves 217 ///@} 218 219 private: 220 221 //---------------------------------------- 222 /**@name Private data */ 223 ///@{ 224 Temp temp; ///< Temporary storage 225 ///@} 226 227 //---------------------------------------- 228 /**@name Solving 229 These helper methods are used during the factorization process. 230 The solve*-methods solve lower and upper triangular systems from 231 the left or from the right, respectively The methods with '2' in 232 the end solve two systems at the same time. The methods with 233 "Eps" in the end consider elements smaller then the passed epsilon 234 as zero. 235 */ 236 ///@{ 237 /// 238 void solveUright(Rational* wrk, Rational* vec); 239 /// 240 int solveUrightEps(Rational* vec, int* nonz, Rational* rhs); 241 /// 242 void solveUright2(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2); 243 /// 244 int solveUright2eps(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2, int* nonz); 245 /// 246 void solveLright2(Rational* vec1, Rational* vec2); 247 /// 248 void solveUpdateRight(Rational* vec); 249 /// 250 void solveUpdateRight2(Rational* vec1, Rational* vec2); 251 /// 252 void solveUleft(Rational* work, Rational* vec); 253 /// 254 void solveUleft2(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2); 255 /// 256 int solveLleft2forest(Rational* vec1, int* /* nonz */, Rational* vec2); 257 /// 258 void solveLleft2(Rational* vec1, int* /* nonz */, Rational* vec2); 259 /// 260 int solveLleftForest(Rational* vec, int* /* nonz */); 261 /// 262 void solveLleft(Rational* vec); 263 /// 264 int solveLleftEps(Rational* vec, int* nonz); 265 /// 266 void solveUpdateLeft(Rational* vec); 267 /// 268 void solveUpdateLeft2(Rational* vec1, Rational* vec2); 269 270 /// 271 int vSolveLright(Rational* vec, int* ridx, int rn); 272 /// 273 void vSolveLright2(Rational* vec, int* ridx, int* rnptr, 274 Rational* vec2, int* ridx2, int* rn2ptr); 275 /// 276 void vSolveLright3(Rational* vec, int* ridx, int* rnptr, 277 Rational* vec2, int* ridx2, int* rn2ptr, 278 Rational* vec3, int* ridx3, int* rn3ptr); 279 /// 280 int vSolveUright(Rational* vec, int* vidx, Rational* rhs, int* ridx, int rn); 281 /// 282 void vSolveUrightNoNZ(Rational* vec, Rational* rhs, int* ridx, int rn); 283 /// 284 int vSolveUright2(Rational* vec, int* vidx, Rational* rhs, int* ridx, int rn, 285 Rational* vec2, Rational* rhs2, int* ridx2, int rn2); 286 /// 287 int vSolveUpdateRight(Rational* vec, int* ridx, int n); 288 /// 289 void vSolveUpdateRightNoNZ(Rational* vec); 290 /// 291 int solveUleft(Rational* vec, int* vecidx, Rational* rhs, int* rhsidx, int rhsn); 292 /// 293 void solveUleftNoNZ(Rational* vec, Rational* rhs, int* rhsidx, int rhsn); 294 /// 295 int solveLleftForest(Rational* vec, int* nonz, int n); 296 /// 297 void solveLleftForestNoNZ(Rational* vec); 298 /// 299 int solveLleft(Rational* vec, int* nonz, int rn); 300 /// 301 void solveLleftNoNZ(Rational* vec); 302 /// 303 int solveUpdateLeft(Rational* vec, int* nonz, int n); 304 305 /// 306 void forestPackColumns(); 307 /// 308 void forestMinColMem(int size); 309 /// 310 void forestReMaxCol(int col, int len); 311 312 /// 313 void initPerm(); 314 /// 315 void initFactorMatrix(const SVectorRational** vec); 316 /// 317 void minLMem(int size); 318 /// 319 void setPivot(const int p_stage, const int p_col, const int p_row, const Rational& val); 320 /// 321 void colSingletons(); 322 /// 323 void rowSingletons(); 324 325 /// 326 void initFactorRings(); 327 /// 328 void freeFactorRings(); 329 330 /// 331 int setupColVals(); 332 /// 333 void setupRowVals(); 334 335 /// 336 void eliminateRowSingletons(); 337 /// 338 void eliminateColSingletons(); 339 /// 340 void selectPivots(const Rational& threshold); 341 /// 342 int updateRow(int r, int lv, int prow, int pcol, const Rational& pval); 343 344 /// 345 void eliminatePivot(int prow, int pos); 346 /// 347 void eliminateNucleus(const Rational& threshold); 348 /// 349 void minRowMem(int size); 350 /// 351 void minColMem(int size); 352 /// 353 void remaxCol(int p_col, int len); 354 /// 355 void packRows(); 356 /// 357 void packColumns(); 358 /// 359 void remaxRow(int p_row, int len); 360 /// 361 int makeLvec(int p_len, int p_row); 362 /// 363 bool timeLimitReached() 364 { 365 if(timeLimit >= 0.0 && factorTime->time() >= timeLimit) 366 { 367 stat = SLinSolverRational::TIME; 368 return true; 369 } 370 else 371 return false; 372 } 373 ///@} 374 375 protected: 376 377 //---------------------------------------- 378 /**@name Solver methods */ 379 ///@{ 380 /// 381 void solveLright(Rational* vec); 382 /// 383 int solveRight4update(Rational* vec, int* nonz, Rational* rhs, 384 Rational* forest, int* forestNum, int* forestIdx); 385 /// 386 void solveRight(Rational* vec, Rational* rhs); 387 /// 388 int solveRight2update(Rational* vec1, Rational* vec2, Rational* rhs1, 389 Rational* rhs2, int* nonz, Rational* forest, int* forestNum, int* forestIdx); 390 /// 391 void solveRight2(Rational* vec1, Rational* vec2, Rational* rhs1, Rational* rhs2); 392 /// 393 void solveLeft(Rational* vec, Rational* rhs); 394 /// 395 int solveLeftEps(Rational* vec, Rational* rhs, int* nonz); 396 /// 397 int solveLeft2(Rational* vec1, int* nonz, Rational* vec2, Rational* rhs1, Rational* rhs2); 398 399 /// 400 int vSolveRight4update( 401 Rational* vec, int* idx, /* result */ 402 Rational* rhs, int* ridx, int rn, /* rhs & Forest */ 403 Rational* forest, int* forestNum, int* forestIdx); 404 /// 405 int vSolveRight4update2( 406 Rational* vec, int* idx, /* result1 */ 407 Rational* rhs, int* ridx, int rn, /* rhs1 */ 408 Rational* vec2, /* result2 */ 409 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ 410 Rational* forest, int* forestNum, int* forestIdx); 411 /// 412 int vSolveRight4update3( 413 Rational* vec, int* idx, /* result1 */ 414 Rational* rhs, int* ridx, int rn, /* rhs1 */ 415 Rational* vec2, /* result2 */ 416 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ 417 Rational* vec3, /* result3 */ 418 Rational* rhs3, int* ridx3, int rn3, /* rhs3 */ 419 Rational* forest, int* forestNum, int* forestIdx); 420 /// 421 void vSolveRightNoNZ(Rational* vec2, /* result2 */ 422 Rational* rhs2, int* ridx2, int rn2); /* rhs2 */ 423 /// 424 int vSolveLeft( 425 Rational* vec, int* idx, /* result */ 426 Rational* rhs, int* ridx, int rn); /* rhs */ 427 /// 428 void vSolveLeftNoNZ( 429 Rational* vec, /* result */ 430 Rational* rhs, int* ridx, int rn); /* rhs */ 431 /// 432 int vSolveLeft2( 433 Rational* vec, int* idx, /* result */ 434 Rational* rhs, int* ridx, int rn, /* rhs */ 435 Rational* vec2, /* result2 */ 436 Rational* rhs2, int* ridx2, int rn2); /* rhs2 */ 437 /// 438 int vSolveLeft3(Rational* vec, int* idx, /* result */ 439 Rational* rhs, int* ridx, int rn, /* rhs */ 440 Rational* vec2, /* result2 */ 441 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ 442 Rational* vec3, /* result3 */ 443 Rational* rhs3, int* ridx3, int rn3); /* rhs3 */ 444 445 void forestUpdate(int col, Rational* work, int num, int* nonz); 446 447 void update(int p_col, Rational* p_work, const int* p_idx, int num); 448 void updateNoClear(int p_col, const Rational* p_work, const int* p_idx, int num); 449 450 /// 451 void factor(const SVectorRational** vec, ///< Array of column vector pointers 452 const Rational& threshold); ///< pivoting threshold 453 ///@} 454 455 //---------------------------------------- 456 /**@name Debugging */ 457 ///@{ 458 /// 459 void dump() const; 460 461 /// 462 bool isConsistent() const; 463 ///@} 464 }; 465 466 } // namespace soplex 467 #include "clufactor_rational.hpp" 468 #endif // _CLUFACTOR_RATIONAL_H_ 469