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
(1) Event rule_of_three_violation: Class "soplex::CLUFactor<double>::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]
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:
(2) Event copy_ctor: User-defined copy constructor.
Also see events: [rule_of_three_violation][copy_assign][remediation]
95   	      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]
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