/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* */ /* This file is part of the class library */ /* SoPlex --- the Sequential object-oriented simPlex. */ /* */ /* Copyright 1996-2022 Zuse Institute Berlin */ /* */ /* Licensed under the Apache License, Version 2.0 (the "License"); */ /* you may not use this file except in compliance with the License. */ /* You may obtain a copy of the License at */ /* */ /* http://www.apache.org/licenses/LICENSE-2.0 */ /* */ /* Unless required by applicable law or agreed to in writing, software */ /* distributed under the License is distributed on an "AS IS" BASIS, */ /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ /* See the License for the specific language governing permissions and */ /* limitations under the License. */ /* */ /* You should have received a copy of the Apache-2.0 license */ /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */ /* */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /**@file clufactor_rational.h * @brief Implementation of sparse LU factorization with Rational precision. */ #ifndef _CLUFACTOR_RATIONAL_H_ #define _CLUFACTOR_RATIONAL_H_ #include "soplex/spxdefines.h" #include "soplex/slinsolver_rational.h" #include "soplex/timer.h" #include "soplex/svector.h" #include "soplex/rational.h" #include "soplex/basevectors.h" #define WITH_L_ROWS 1 namespace soplex { /**@brief Implementation of sparse LU factorization with Rational precision. * @ingroup Algo * * This class implements a sparse LU factorization with either * FOREST-TOMLIN or ETA updates, using dynamic Markowitz pivoting. */ class CLUFactorRational { public: //---------------------------------------- /**@name Public types */ ///@{ /** Doubly linked ring structure for garbage collection of column or * row file in working matrix. */ struct Dring { Dring* next; Dring* prev; int idx; }; /// Pivot Ring class Pring { public: Pring* next; ///< Pring* prev; ///< int idx; ///< index of pivot row int pos; ///< position of pivot column in row int mkwtz; ///< markowitz number of pivot Pring() ///< constructor : next(0) , prev(0) , idx(0) , pos(0) , mkwtz(0) {} private: Pring(const Pring&); ///< blocked copy constructor Pring& operator= (const Pring&); ///< blocked assignment operator }; ///@} protected: //---------------------------------------- /**@name Protected types */ ///@{ /// Temporary data structures. class Temp { public: int* s_mark; ///< marker VectorRational s_max; ///< maximum absolute value per row (or -1) int* s_cact; ///< lengths of columns of active submatrix int stage; ///< stage of the structure Pring pivots; ///< ring of selected pivot rows Pring* pivot_col; ///< column index handlers for Real linked list Pring* pivot_colNZ; ///< lists for columns to number of nonzeros Pring* pivot_row; ///< row index handlers for Real linked list Pring* pivot_rowNZ; ///< lists for rows to number of nonzeros Temp(); ///< constructor ~Temp(); ///< destructor void init(int p_dim); ///< initialization void clear(); ///< clears the structure private: Temp(const Temp&); ///< blocked copy constructor Temp& operator= (const Temp&); ///< blocked assignment operator }; /// Data structures for saving the row and column permutations. struct Perm { int* orig; ///< orig[p] original index from p int* perm; ///< perm[i] permuted index from i }; /// Data structures for saving the working matrix and U factor. struct U { /// struct Row { Dring list; /*!< \brief Double linked ringlist of vector indices in the order they appear in the row file */ Dring* elem; ///< %Array of ring elements. int used; ///< used entries of arrays idx and val VectorRational val; ///< hold nonzero values int* idx; ///< array of length val.dim() to hold column indices of nonzeros in val int* start; ///< starting positions in val and idx int* len; ///< used nonzeros per row vectors int* max; /*!< \brief maximum available nonzeros per row: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]. */ } row; /// struct Col { Dring list; /*!< \brief Double linked ringlist of vector indices in the order they appear in the column file */ Dring* elem; ///< %Array of ring elements. int size; ///< size of array idx int used; ///< used entries of array idx int* idx; ///< hold row indices of nonzeros VectorRational val; /*!< \brief hold nonzero values: this is only initialized in the end of the factorization with DEFAULT updates. */ int* start; ///< starting positions in val and idx int* len; ///< used nonzeros per column vector int* max; /*!< \brief maximum available nonzeros per colunn: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]. */ } col; }; /// Data structures for saving the working matrix and L factor. struct L { VectorRational val; ///< values of L vectors int* idx; ///< array of size val.dim() storing indices of L vectors int startSize; ///< size of array start int firstUpdate; ///< number of first update L vector int firstUnused; ///< number of first unused L vector int* start; ///< starting positions in val and idx int* row; ///< column indices of L vectors int updateType; ///< type of updates to be used. /* The following arrays have length |firstUpdate|, since they keep * rows of the L-vectors occuring during the factorization (without * updates), only: */ VectorRational rval; ///< values of rows of L int* ridx; ///< indices of rows of L int* rbeg; ///< start of rows in rval and ridx int* rorig; ///< original row permutation int* rperm; ///< original row permutation }; ///@} //---------------------------------------- /**@name Protected data */ ///@{ SLinSolverRational::Status stat; ///< Status indicator. int thedim; ///< dimension of factorized matrix int nzCnt; ///< number of nonzeros in U Rational initMaxabs; ///< maximum abs number in initail Matrix Rational maxabs; ///< maximum abs number in L and U Real rowMemMult; ///< factor of minimum Memory * number of nonzeros Real colMemMult; ///< factor of minimum Memory * number of nonzeros Real lMemMult; ///< factor of minimum Memory * number of nonzeros Perm row; ///< row permutation matrices Perm col; ///< column permutation matrices L l; ///< L matrix VectorRational diag; ///< Array of pivot elements U u; ///< U matrix Rational* work; ///< Working array: must always be left as 0! Timer* factorTime; ///< Time spent in factorizations int factorCount; ///< Number of factorizations Real timeLimit; ///< Time limit on factorization or solves ///@} private: //---------------------------------------- /**@name Private data */ ///@{ Temp temp; ///< Temporary storage ///@} //---------------------------------------- /**@name Solving These helper methods are used during the factorization process. The solve*-methods solve lower and upper triangular systems from the left or from the right, respectively The methods with '2' in the end solve two systems at the same time. The methods with "Eps" in the end consider elements smaller then the passed epsilon as zero. */ ///@{ /// void solveUright(Rational* wrk, Rational* vec); /// int solveUrightEps(Rational* vec, int* nonz, Rational* rhs); /// void solveUright2(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2); /// int solveUright2eps(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2, int* nonz); /// void solveLright2(Rational* vec1, Rational* vec2); /// void solveUpdateRight(Rational* vec); /// void solveUpdateRight2(Rational* vec1, Rational* vec2); /// void solveUleft(Rational* work, Rational* vec); /// void solveUleft2(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2); /// int solveLleft2forest(Rational* vec1, int* /* nonz */, Rational* vec2); /// void solveLleft2(Rational* vec1, int* /* nonz */, Rational* vec2); /// int solveLleftForest(Rational* vec, int* /* nonz */); /// void solveLleft(Rational* vec); /// int solveLleftEps(Rational* vec, int* nonz); /// void solveUpdateLeft(Rational* vec); /// void solveUpdateLeft2(Rational* vec1, Rational* vec2); /// int vSolveLright(Rational* vec, int* ridx, int rn); /// void vSolveLright2(Rational* vec, int* ridx, int* rnptr, Rational* vec2, int* ridx2, int* rn2ptr); /// void vSolveLright3(Rational* vec, int* ridx, int* rnptr, Rational* vec2, int* ridx2, int* rn2ptr, Rational* vec3, int* ridx3, int* rn3ptr); /// int vSolveUright(Rational* vec, int* vidx, Rational* rhs, int* ridx, int rn); /// void vSolveUrightNoNZ(Rational* vec, Rational* rhs, int* ridx, int rn); /// int vSolveUright2(Rational* vec, int* vidx, Rational* rhs, int* ridx, int rn, Rational* vec2, Rational* rhs2, int* ridx2, int rn2); /// int vSolveUpdateRight(Rational* vec, int* ridx, int n); /// void vSolveUpdateRightNoNZ(Rational* vec); /// int solveUleft(Rational* vec, int* vecidx, Rational* rhs, int* rhsidx, int rhsn); /// void solveUleftNoNZ(Rational* vec, Rational* rhs, int* rhsidx, int rhsn); /// int solveLleftForest(Rational* vec, int* nonz, int n); /// void solveLleftForestNoNZ(Rational* vec); /// int solveLleft(Rational* vec, int* nonz, int rn); /// void solveLleftNoNZ(Rational* vec); /// int solveUpdateLeft(Rational* vec, int* nonz, int n); /// void forestPackColumns(); /// void forestMinColMem(int size); /// void forestReMaxCol(int col, int len); /// void initPerm(); /// void initFactorMatrix(const SVectorRational** vec); /// void minLMem(int size); /// void setPivot(const int p_stage, const int p_col, const int p_row, const Rational& val); /// void colSingletons(); /// void rowSingletons(); /// void initFactorRings(); /// void freeFactorRings(); /// int setupColVals(); /// void setupRowVals(); /// void eliminateRowSingletons(); /// void eliminateColSingletons(); /// void selectPivots(const Rational& threshold); /// int updateRow(int r, int lv, int prow, int pcol, const Rational& pval); /// void eliminatePivot(int prow, int pos); /// void eliminateNucleus(const Rational& threshold); /// void minRowMem(int size); /// void minColMem(int size); /// void remaxCol(int p_col, int len); /// void packRows(); /// void packColumns(); /// void remaxRow(int p_row, int len); /// int makeLvec(int p_len, int p_row); /// bool timeLimitReached() { if(timeLimit >= 0.0 && factorTime->time() >= timeLimit) { stat = SLinSolverRational::TIME; return true; } else return false; } ///@} protected: //---------------------------------------- /**@name Solver methods */ ///@{ /// void solveLright(Rational* vec); /// int solveRight4update(Rational* vec, int* nonz, Rational* rhs, Rational* forest, int* forestNum, int* forestIdx); /// void solveRight(Rational* vec, Rational* rhs); /// int solveRight2update(Rational* vec1, Rational* vec2, Rational* rhs1, Rational* rhs2, int* nonz, Rational* forest, int* forestNum, int* forestIdx); /// void solveRight2(Rational* vec1, Rational* vec2, Rational* rhs1, Rational* rhs2); /// void solveLeft(Rational* vec, Rational* rhs); /// int solveLeftEps(Rational* vec, Rational* rhs, int* nonz); /// int solveLeft2(Rational* vec1, int* nonz, Rational* vec2, Rational* rhs1, Rational* rhs2); /// int vSolveRight4update( Rational* vec, int* idx, /* result */ Rational* rhs, int* ridx, int rn, /* rhs & Forest */ Rational* forest, int* forestNum, int* forestIdx); /// int vSolveRight4update2( Rational* vec, int* idx, /* result1 */ Rational* rhs, int* ridx, int rn, /* rhs1 */ Rational* vec2, /* result2 */ Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ Rational* forest, int* forestNum, int* forestIdx); /// int vSolveRight4update3( Rational* vec, int* idx, /* result1 */ Rational* rhs, int* ridx, int rn, /* rhs1 */ Rational* vec2, /* result2 */ Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ Rational* vec3, /* result3 */ Rational* rhs3, int* ridx3, int rn3, /* rhs3 */ Rational* forest, int* forestNum, int* forestIdx); /// void vSolveRightNoNZ(Rational* vec2, /* result2 */ Rational* rhs2, int* ridx2, int rn2); /* rhs2 */ /// int vSolveLeft( Rational* vec, int* idx, /* result */ Rational* rhs, int* ridx, int rn); /* rhs */ /// void vSolveLeftNoNZ( Rational* vec, /* result */ Rational* rhs, int* ridx, int rn); /* rhs */ /// int vSolveLeft2( Rational* vec, int* idx, /* result */ Rational* rhs, int* ridx, int rn, /* rhs */ Rational* vec2, /* result2 */ Rational* rhs2, int* ridx2, int rn2); /* rhs2 */ /// int vSolveLeft3(Rational* vec, int* idx, /* result */ Rational* rhs, int* ridx, int rn, /* rhs */ Rational* vec2, /* result2 */ Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ Rational* vec3, /* result3 */ Rational* rhs3, int* ridx3, int rn3); /* rhs3 */ void forestUpdate(int col, Rational* work, int num, int* nonz); void update(int p_col, Rational* p_work, const int* p_idx, int num); void updateNoClear(int p_col, const Rational* p_work, const int* p_idx, int num); /// void factor(const SVectorRational** vec, ///< Array of column vector pointers const Rational& threshold); ///< pivoting threshold ///@} //---------------------------------------- /**@name Debugging */ ///@{ /// void dump() const; /// bool isConsistent() const; ///@} }; } // namespace soplex #include "clufactor_rational.hpp" #endif // _CLUFACTOR_RATIONAL_H_