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
(1) Event rule_of_three_violation: Class "soplex::CLUFactorRational::Pring" has a user definition for at least one special function (copy constructor, copy assignment, destructor) but not all. If one of these functions requires a user definition then the others likely do as well.
(4) Event remediation: Add user-definition for a destructor.
Also see events: [copy_ctor][copy_assign]
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:
(2) Event copy_ctor: User-defined copy constructor.
Also see events: [rule_of_three_violation][copy_assign][remediation]
84   	      Pring(const Pring&);             ///< blocked copy constructor
(3) Event copy_assign: User-defined copy assignment operator.
Also see events: [rule_of_three_violation][copy_ctor][remediation]
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