1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_sparse_decomposition_h
17 #define dealii_sparse_decomposition_h
18
19 #include <deal.II/base/config.h>
20
21 #include <deal.II/lac/sparse_matrix.h>
22
23 #include <cmath>
24
25 DEAL_II_NAMESPACE_OPEN
26
27 /*! @addtogroup Preconditioners
28 *@{
29 */
30
31 /**
32 * Abstract base class for incomplete decompositions of a sparse matrix into
33 * sparse factors. This class can't be used by itself, but only as the base
34 * class of derived classes that actually implement particular decompositions
35 * such as SparseILU or SparseMIC.
36 *
37 * The decomposition is stored as a sparse matrix which is why this class is
38 * derived from the SparseMatrix. Since it is not a matrix in the usual sense
39 * (the stored entries are not those of a matrix, but of the two factors of
40 * the original matrix), the derivation is <tt>protected</tt> rather than
41 * <tt>public</tt>.
42 *
43 *
44 * <h3>Fill-in</h3>
45 *
46 * Sparse decompositions are frequently used with additional fill-in, i.e.,
47 * the sparsity structure of the decomposition is denser than that of the
48 * matrix to be decomposed. The initialize() function of this class allows
49 * this fill-in via the AdditionalData object as long as all entries present
50 * in the original matrix are present in the decomposition also, i.e. the
51 * sparsity pattern of the decomposition is a superset of the sparsity pattern
52 * in the original matrix.
53 *
54 * Such fill-in can be accomplished by various ways, one of which is the copy-
55 * constructor of the SparsityPattern class that allows the addition of side-
56 * diagonals to a given sparsity structure.
57 *
58 *
59 * <h3>Unified use of preconditioners</h3>
60 *
61 * While objects of this class can not be used directly (this class is only a
62 * base class for others implementing actual decompositions), derived classes
63 * such as SparseILU and SparseMIC can be used in the usual form as
64 * preconditioners. For example, this works:
65 * @code
66 * SparseILU<double> ilu;
67 * ilu.initialize(matrix, SparseILU<double>::AdditionalData(...));
68 *
69 * somesolver.solve (A, x, f, ilu);
70 * @endcode
71 *
72 * Through the AdditionalData object it is possible to specify additional
73 * parameters of the LU decomposition.
74 *
75 * 1/ The matrix diagonal can be strengthened by adding
76 * <code>strengthen_diagonal</code> times the sum of the absolute row entries
77 * of each row to the respective diagonal entries. By default no strengthening
78 * is performed.
79 *
80 * 2/ By default, each initialize() function call creates its own sparsity.
81 * For that, it copies the sparsity of <code>matrix</code> and adds a specific
82 * number of extra off diagonal entries specified by
83 * <code>extra_off_diagonals</code>.
84 *
85 * 3/ By setting <code>use_previous_sparsity=true</code> the sparsity is not
86 * recreated but the sparsity of the previous initialize() call is reused
87 * (recycled). This might be useful when several linear problems on the same
88 * sparsity need to solved, as for example several Newton iteration steps on
89 * the same triangulation. The default is <code>false</code>.
90 *
91 * 4/ It is possible to give a user defined sparsity to
92 * <code>use_this_sparsity</code>. Then, no sparsity is created but
93 * <code>*use_this_sparsity</code> is used to store the decomposed matrix. For
94 * restrictions on the sparsity see section `Fill-in' above).
95 *
96 *
97 * <h3>Particular implementations</h3>
98 *
99 * It is enough to override the initialize() and vmult() methods to implement
100 * particular LU decompositions, like the true LU, or the Cholesky
101 * decomposition. Additionally, if that decomposition needs fine tuned
102 * diagonal strengthening on a per row basis, it may override the
103 * get_strengthen_diagonal() method.
104 */
105 template <typename number>
106 class SparseLUDecomposition : protected SparseMatrix<number>,
107 public virtual Subscriptor
108 {
109 protected:
110 /**
111 * Constructor. Does nothing.
112 *
113 * Call the initialize() function before using this object as preconditioner
114 * (vmult()).
115 */
116 SparseLUDecomposition();
117
118 public:
119 /**
120 * Declare type for container size.
121 */
122 using size_type = typename SparseMatrix<number>::size_type;
123
124 /**
125 * Destruction. Mark the destructor pure to ensure that this class isn't
126 * used directly, but only its derived classes.
127 */
128 virtual ~SparseLUDecomposition() override = 0;
129
130 /**
131 * Deletes all member variables. Leaves the class in the state that it had
132 * directly after calling the constructor
133 */
134 virtual void
135 clear() override;
136
137 /**
138 * Parameters for SparseDecomposition.
139 */
140 class AdditionalData
141 {
142 public:
143 /**
144 * Constructor. For the parameters' description, see below.
145 */
146 AdditionalData(const double strengthen_diagonal = 0,
147 const unsigned int extra_off_diagonals = 0,
148 const bool use_previous_sparsity = false,
149 const SparsityPattern *use_this_sparsity = nullptr);
150
151 /**
152 * <code>strengthen_diag</code> times the sum of absolute row entries is
153 * added to the diagonal entries.
154 *
155 * Per default, this value is zero, i.e. the diagonal is not strengthened.
156 */
157 double strengthen_diagonal;
158
159 /**
160 * By default, the <code>initialize(matrix, data)</code> function creates
161 * its own sparsity. This sparsity has the same SparsityPattern as
162 * <code>matrix</code> with some extra off diagonals the number of which
163 * is specified by <code>extra_off_diagonals</code>.
164 *
165 * The user can give a SparsityPattern to <code>use_this_sparsity</code>.
166 * Then this sparsity is used and the <code>extra_off_diagonals</code>
167 * argument is ignored.
168 */
169 unsigned int extra_off_diagonals;
170
171 /**
172 * If this flag is true the initialize() function uses the same sparsity
173 * that was used during the previous initialize() call.
174 *
175 * This might be useful when several linear problems on the same sparsity
176 * need to solved, as for example several Newton iteration steps on the
177 * same triangulation.
178 */
179 bool use_previous_sparsity;
180
181 /**
182 * When a SparsityPattern is given to this argument, the initialize()
183 * function calls <code>reinit(*use_this_sparsity)</code> causing this
184 * sparsity to be used.
185 *
186 * Note that the sparsity structures of <code>*use_this_sparsity</code>
187 * and the matrix passed to the initialize function need not be equal.
188 * Fill-in is allowed, as well as filtering out some elements in the
189 * matrix.
190 */
191 const SparsityPattern *use_this_sparsity;
192 };
193
194 /**
195 * This function needs to be called before an object of this class is used
196 * as preconditioner.
197 *
198 * For more detail about possible parameters, see the class documentation
199 * and the documentation of the SparseLUDecomposition::AdditionalData class.
200 *
201 * According to the <code>parameters</code>, this function creates a new
202 * SparsityPattern or keeps the previous sparsity or takes the sparsity
203 * given by the user to <code>data</code>. Then, this function performs the
204 * LU decomposition.
205 *
206 * After this function is called the preconditioner is ready to be used
207 * (using the <code>vmult</code> function of derived classes).
208 */
209 template <typename somenumber>
210 void
211 initialize(const SparseMatrix<somenumber> &matrix,
212 const AdditionalData parameters);
213
214 /**
215 * Return whether the object is empty. It calls the inherited
216 * SparseMatrix::empty() function.
217 */
218 bool
219 empty() const;
220
221 /**
222 * Return the dimension of the codomain (or range) space. It calls the
223 * inherited SparseMatrix::m() function. Note that the matrix is of
224 * dimension $m \times n$.
225 */
226 size_type
227 m() const;
228
229 /**
230 * Return the dimension of the domain space. It calls the inherited
231 * SparseMatrix::n() function. Note that the matrix is of dimension $m
232 * \times n$.
233 */
234 size_type
235 n() const;
236
237 /**
238 * Adding Matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i> with
239 * <i>M</i> being this matrix.
240 *
241 * Source and destination must not be the same vector.
242 */
243 template <class OutVector, class InVector>
244 void
245 vmult_add(OutVector &dst, const InVector &src) const;
246
247 /**
248 * Adding Matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
249 * <i>dst</i> with <i>M</i> being this matrix. This function does the same
250 * as vmult_add() but takes the transposed matrix.
251 *
252 * Source and destination must not be the same vector.
253 */
254 template <class OutVector, class InVector>
255 void
256 Tvmult_add(OutVector &dst, const InVector &src) const;
257
258 /**
259 * Determine an estimate for the memory consumption (in bytes) of this
260 * object.
261 */
262 virtual std::size_t
263 memory_consumption() const;
264
265 /**
266 * @addtogroup Exceptions
267 * @{
268 */
269
270 /**
271 * Exception
272 */
273 DeclException1(ExcInvalidStrengthening,
274 double,
275 << "The strengthening parameter " << arg1
276 << " is not greater or equal than zero!");
277 //@}
278 protected:
279 /**
280 * Copies the passed SparseMatrix onto this object. This object's sparsity
281 * pattern remains unchanged.
282 */
283 template <typename somenumber>
284 void
285 copy_from(const SparseMatrix<somenumber> &matrix);
286
287 /**
288 * Performs the strengthening loop. For each row calculates the sum of
289 * absolute values of its elements, determines the strengthening factor
290 * (through get_strengthen_diagonal()) sf and multiplies the diagonal entry
291 * with <code>sf+1</code>.
292 */
293 virtual void
294 strengthen_diagonal_impl();
295
296 /**
297 * In the decomposition phase, computes a strengthening factor for the
298 * diagonal entry in row <code>row</code> with sum of absolute values of its
299 * elements <code>rowsum</code>.
300 *
301 * @note The default implementation in SparseLUDecomposition returns
302 * <code>strengthen_diagonal</code>'s value. This variable is set to
303 * a nonzero value in several of the derived classes.
304 */
305 virtual number
306 get_strengthen_diagonal(const number rowsum, const size_type row) const;
307
308 /**
309 * The default strengthening value, returned by get_strengthen_diagonal().
310 */
311 double strengthen_diagonal;
312
313 /**
314 * For every row in the underlying SparsityPattern, this array contains a
315 * pointer to the row's first afterdiagonal entry. Becomes available after
316 * invocation of prebuild_lower_bound().
317 */
318 std::vector<const size_type *> prebuilt_lower_bound;
319
320 /**
321 * Fills the #prebuilt_lower_bound array.
322 */
323 void
324 prebuild_lower_bound();
325
326 private:
327 /**
328 * In general this pointer is zero except for the case that no
329 * SparsityPattern is given to this class. Then, a SparsityPattern is
330 * created and is passed down to the SparseMatrix base class.
331 *
332 * Nevertheless, the SparseLUDecomposition needs to keep ownership of this
333 * sparsity. It keeps this pointer to it enabling it to delete this sparsity
334 * at destruction time.
335 */
336 SparsityPattern *own_sparsity;
337 };
338
339 /*@}*/
340 //---------------------------------------------------------------------------
341
342 #ifndef DOXYGEN
343
344 template <typename number>
345 inline number
get_strengthen_diagonal(const number,const size_type)346 SparseLUDecomposition<number>::get_strengthen_diagonal(
347 const number /*rowsum*/,
348 const size_type /*row*/) const
349 {
350 return strengthen_diagonal;
351 }
352
353
354
355 template <typename number>
356 inline bool
empty()357 SparseLUDecomposition<number>::empty() const
358 {
359 return SparseMatrix<number>::empty();
360 }
361
362
363 template <typename number>
364 inline typename SparseLUDecomposition<number>::size_type
m()365 SparseLUDecomposition<number>::m() const
366 {
367 return SparseMatrix<number>::m();
368 }
369
370
371 template <typename number>
372 inline typename SparseLUDecomposition<number>::size_type
n()373 SparseLUDecomposition<number>::n() const
374 {
375 return SparseMatrix<number>::n();
376 }
377
378 // Note: This function is required for full compatibility with
379 // the LinearOperator class. ::MatrixInterfaceWithVmultAdd
380 // picks up the vmult_add function in the protected SparseMatrix
381 // base class.
382 template <typename number>
383 template <class OutVector, class InVector>
384 inline void
vmult_add(OutVector & dst,const InVector & src)385 SparseLUDecomposition<number>::vmult_add(OutVector & dst,
386 const InVector &src) const
387 {
388 OutVector tmp;
389 tmp.reinit(dst);
390 this->vmult(tmp, src);
391 dst += tmp;
392 }
393
394 // Note: This function is required for full compatibility with
395 // the LinearOperator class. ::MatrixInterfaceWithVmultAdd
396 // picks up the vmult_add function in the protected SparseMatrix
397 // base class.
398 template <typename number>
399 template <class OutVector, class InVector>
400 inline void
Tvmult_add(OutVector & dst,const InVector & src)401 SparseLUDecomposition<number>::Tvmult_add(OutVector & dst,
402 const InVector &src) const
403 {
404 OutVector tmp;
405 tmp.reinit(dst);
406 this->Tvmult(tmp, src);
407 dst += tmp;
408 }
409
410 //---------------------------------------------------------------------------
411
412
413 template <typename number>
AdditionalData(const double strengthen_diag,const unsigned int extra_off_diag,const bool use_prev_sparsity,const SparsityPattern * use_this_spars)414 SparseLUDecomposition<number>::AdditionalData::AdditionalData(
415 const double strengthen_diag,
416 const unsigned int extra_off_diag,
417 const bool use_prev_sparsity,
418 const SparsityPattern *use_this_spars)
419 : strengthen_diagonal(strengthen_diag)
420 , extra_off_diagonals(extra_off_diag)
421 , use_previous_sparsity(use_prev_sparsity)
422 , use_this_sparsity(use_this_spars)
423 {}
424
425
426 #endif // DOXYGEN
427
428 DEAL_II_NAMESPACE_CLOSE
429
430 #endif // dealii_sparse_decomposition_h
431