1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 1999 - 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 17 #ifndef dealii_sparse_ilu_h 18 #define dealii_sparse_ilu_h 19 20 21 #include <deal.II/base/config.h> 22 23 #include <deal.II/lac/exceptions.h> 24 #include <deal.II/lac/sparse_decomposition.h> 25 #include <deal.II/lac/sparse_matrix.h> 26 27 DEAL_II_NAMESPACE_OPEN 28 29 /*! @addtogroup Preconditioners 30 *@{ 31 */ 32 33 /** 34 * This class computes an Incomplete LU (ILU) decomposition of a sparse 35 * matrix, using either the same sparsity pattern or a different one. By 36 * incomplete we mean that unlike the exact decomposition, the incomplete one 37 * is also computed using sparse factors, and entries in the decomposition 38 * that do not fit into the given sparsity structure are discarded. 39 * 40 * The algorithm used by this class is essentially a copy of the algorithm 41 * given in the book Y. Saad: "Iterative methods for sparse linear systems", 42 * second edition, in section 10.3.2. 43 * 44 * 45 * <h3>Usage and state management</h3> 46 * 47 * Refer to SparseLUDecomposition documentation for suggested usage and state 48 * management. This class is used in the 49 * @ref step_22 "step-22" 50 * tutorial program. 51 * 52 * @note Instantiations for this template are provided for <tt>@<float@> and 53 * @<double@></tt>; others can be generated in application programs (see the 54 * section on 55 * @ref Instantiations 56 * in the manual). 57 */ 58 template <typename number> 59 class SparseILU : public SparseLUDecomposition<number> 60 { 61 public: 62 /** 63 * Declare type for container size. 64 */ 65 using size_type = typename SparseLUDecomposition<number>::size_type; 66 67 /** 68 * Constructor. Does nothing. 69 * 70 * Call the @p initialize function before using this object as 71 * preconditioner. 72 */ 73 SparseILU() = default; 74 75 /** 76 * Make SparseLUDecomposition::AdditionalData accessible to this class as 77 * well. 78 */ 79 using AdditionalData = typename SparseLUDecomposition<number>::AdditionalData; 80 81 /** 82 * Perform the incomplete LU factorization of the given matrix. 83 * 84 * This function needs to be called before an object of this class is used 85 * as preconditioner. 86 * 87 * For more details about possible parameters, see the class documentation 88 * of SparseLUDecomposition and the documentation of the @p 89 * SparseLUDecomposition::AdditionalData class. 90 * 91 * According to the @p parameters, this function creates a new 92 * SparsityPattern or keeps the previous sparsity or takes the sparsity 93 * given by the user to @p data. Then, this function performs the LU 94 * decomposition. 95 * 96 * After this function is called the preconditioner is ready to be used. 97 */ 98 template <typename somenumber> 99 void 100 initialize(const SparseMatrix<somenumber> &matrix, 101 const AdditionalData & parameters = AdditionalData()); 102 103 /** 104 * Apply the incomplete decomposition, i.e. do one forward-backward step 105 * $dst=(LU)^{-1}src$. 106 * 107 * The initialize() function needs to be called before. 108 */ 109 template <typename somenumber> 110 void 111 vmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const; 112 113 114 /** 115 * Apply the transpose of the incomplete decomposition, i.e. do one forward- 116 * backward step $dst=(LU)^{-T}src$. 117 * 118 * The initialize() function needs to be called before. 119 */ 120 template <typename somenumber> 121 void 122 Tvmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const; 123 124 125 /** 126 * Determine an estimate for the memory consumption (in bytes) of this 127 * object. 128 */ 129 std::size_t 130 memory_consumption() const override; 131 132 /** 133 * @addtogroup Exceptions 134 * @{ 135 */ 136 137 /** 138 * Exception 139 */ 140 DeclException1(ExcInvalidStrengthening, 141 double, 142 << "The strengthening parameter " << arg1 143 << " is not greater or equal than zero!"); 144 /** 145 * Exception 146 */ 147 DeclException1(ExcZeroPivot, 148 size_type, 149 << "While computing the ILU decomposition, the algorithm " 150 "found a zero pivot on the diagonal of row " 151 << arg1 152 << ". This must stop the ILU algorithm because it means " 153 "that the matrix for which you try to compute a " 154 "decomposition is singular."); 155 //@} 156 }; 157 158 /*@}*/ 159 //--------------------------------------------------------------------------- 160 161 162 DEAL_II_NAMESPACE_CLOSE 163 164 #endif // dealii_sparse_ilu_h 165