1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2004 - 2020 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_petsc_block_sparse_matrix_h 17 #define dealii_petsc_block_sparse_matrix_h 18 19 20 #include <deal.II/base/config.h> 21 22 #ifdef DEAL_II_WITH_PETSC 23 24 # include <deal.II/lac/block_matrix_base.h> 25 # include <deal.II/lac/block_sparsity_pattern.h> 26 # include <deal.II/lac/exceptions.h> 27 # include <deal.II/lac/petsc_block_vector.h> 28 # include <deal.II/lac/petsc_sparse_matrix.h> 29 30 # include <cmath> 31 32 DEAL_II_NAMESPACE_OPEN 33 34 35 36 namespace PETScWrappers 37 { 38 namespace MPI 39 { 40 /*! @addtogroup PETScWrappers 41 *@{ 42 */ 43 44 /** 45 * Blocked sparse matrix based on the PETScWrappers::MPI::SparseMatrix 46 * class. This class implements the functions that are specific to the 47 * PETSc SparseMatrix base objects for a blocked sparse matrix, and leaves 48 * the actual work relaying most of the calls to the individual blocks to 49 * the functions implemented in the base class. See there also for a 50 * description of when this class is useful. 51 * 52 * In contrast to the deal.II-type SparseMatrix class, the PETSc matrices 53 * do not have external objects for the sparsity patterns. Thus, one does 54 * not determine the size of the individual blocks of a block matrix of 55 * this type by attaching a block sparsity pattern, but by calling 56 * reinit() to set the number of blocks and then by setting the size of 57 * each block separately. In order to fix the data structures of the block 58 * matrix, it is then necessary to let it know that we have changed the 59 * sizes of the underlying matrices. For this, one has to call the 60 * collect_sizes() function, for much the same reason as is documented 61 * with the BlockSparsityPattern class. 62 * 63 * @ingroup Matrix1 @see 64 * @ref GlossBlockLA "Block (linear algebra)" 65 */ 66 class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix> 67 { 68 public: 69 /** 70 * Typedef the base class for simpler access to its own alias. 71 */ 72 using BaseClass = BlockMatrixBase<SparseMatrix>; 73 74 /** 75 * Typedef the type of the underlying matrix. 76 */ 77 using BlockType = BaseClass::BlockType; 78 79 /** 80 * Import the alias from the base class. 81 */ 82 using value_type = BaseClass::value_type; 83 using pointer = BaseClass::pointer; 84 using const_pointer = BaseClass::const_pointer; 85 using reference = BaseClass::reference; 86 using const_reference = BaseClass::const_reference; 87 using size_type = BaseClass::size_type; 88 using iterator = BaseClass::iterator; 89 using const_iterator = BaseClass::const_iterator; 90 91 /** 92 * Constructor; initializes the matrix to be empty, without any 93 * structure, i.e. the matrix is not usable at all. This constructor is 94 * therefore only useful for matrices which are members of a class. All 95 * other matrices should be created at a point in the data flow where 96 * all necessary information is available. 97 * 98 * You have to initialize the matrix before usage with 99 * reinit(BlockSparsityPattern). The number of blocks per row and column 100 * are then determined by that function. 101 */ 102 BlockSparseMatrix() = default; 103 104 /** 105 * Destructor. 106 */ 107 ~BlockSparseMatrix() override = default; 108 109 /** 110 * Pseudo copy operator only copying empty objects. The sizes of the 111 * block matrices need to be the same. 112 */ 113 BlockSparseMatrix & 114 operator=(const BlockSparseMatrix &); 115 116 /** 117 * This operator assigns a scalar to a matrix. Since this does usually 118 * not make much sense (should we set all matrix entries to this value? 119 * Only the nonzero entries of the sparsity pattern?), this operation is 120 * only allowed if the actual value to be assigned is zero. This 121 * operator only exists to allow for the obvious notation 122 * <tt>matrix=0</tt>, which sets all elements of the matrix to zero, but 123 * keep the sparsity pattern previously used. 124 */ 125 BlockSparseMatrix & 126 operator=(const double d); 127 128 /** 129 * Resize the matrix, by setting the number of block rows and columns. 130 * This deletes all blocks and replaces them with uninitialized ones, 131 * i.e. ones for which also the sizes are not yet set. You have to do 132 * that by calling the @p reinit functions of the blocks themselves. Do 133 * not forget to call collect_sizes() after that on this object. 134 * 135 * The reason that you have to set sizes of the blocks yourself is that 136 * the sizes may be varying, the maximum number of elements per row may 137 * be varying, etc. It is simpler not to reproduce the interface of the 138 * SparsityPattern class here but rather let the user call whatever 139 * function they desire. 140 */ 141 void 142 reinit(const size_type n_block_rows, const size_type n_block_columns); 143 144 145 /** 146 * Efficiently reinit the block matrix for a parallel computation. Only 147 * the BlockSparsityPattern of the Simple type can efficiently store 148 * large sparsity patterns in parallel, so this is the only supported 149 * argument. The IndexSets describe the locally owned range of DoFs for 150 * each block. Note that the IndexSets needs to be ascending and 1:1. 151 * For a symmetric structure hand in the same vector for the first two 152 * arguments. 153 */ 154 void 155 reinit(const std::vector<IndexSet> & rows, 156 const std::vector<IndexSet> & cols, 157 const BlockDynamicSparsityPattern &bdsp, 158 const MPI_Comm & com); 159 160 161 /** 162 * Same as above but for a symmetric structure only. 163 */ 164 void 165 reinit(const std::vector<IndexSet> & sizes, 166 const BlockDynamicSparsityPattern &bdsp, 167 const MPI_Comm & com); 168 169 170 171 /** 172 * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this 173 * matrix. 174 */ 175 void 176 vmult(BlockVector &dst, const BlockVector &src) const; 177 178 /** 179 * Matrix-vector multiplication. Just like the previous function, but 180 * only applicable if the matrix has only one block column. 181 */ 182 void 183 vmult(BlockVector &dst, const Vector &src) const; 184 185 /** 186 * Matrix-vector multiplication. Just like the previous function, but 187 * only applicable if the matrix has only one block row. 188 */ 189 void 190 vmult(Vector &dst, const BlockVector &src) const; 191 192 /** 193 * Matrix-vector multiplication. Just like the previous function, but 194 * only applicable if the matrix has only one block. 195 */ 196 void 197 vmult(Vector &dst, const Vector &src) const; 198 199 /** 200 * Matrix-vector multiplication: let $dst = M^T*src$ with $M$ being this 201 * matrix. This function does the same as vmult() but takes the 202 * transposed matrix. 203 */ 204 void 205 Tvmult(BlockVector &dst, const BlockVector &src) const; 206 207 /** 208 * Matrix-vector multiplication. Just like the previous function, but 209 * only applicable if the matrix has only one block row. 210 */ 211 void 212 Tvmult(BlockVector &dst, const Vector &src) const; 213 214 /** 215 * Matrix-vector multiplication. Just like the previous function, but 216 * only applicable if the matrix has only one block column. 217 */ 218 void 219 Tvmult(Vector &dst, const BlockVector &src) const; 220 221 /** 222 * Matrix-vector multiplication. Just like the previous function, but 223 * only applicable if the matrix has only one block. 224 */ 225 void 226 Tvmult(Vector &dst, const Vector &src) const; 227 228 /** 229 * This function collects the sizes of the sub-objects and stores them 230 * in internal arrays, in order to be able to relay global indices into 231 * the matrix to indices into the subobjects. You *must* call this 232 * function each time after you have changed the size of the sub- 233 * objects. 234 */ 235 void 236 collect_sizes(); 237 238 /** 239 * Return the partitioning of the domain space of this matrix, i.e., the 240 * partitioning of the vectors this matrix has to be multiplied with. 241 */ 242 std::vector<IndexSet> 243 locally_owned_domain_indices() const; 244 245 /** 246 * Return the partitioning of the range space of this matrix, i.e., the 247 * partitioning of the vectors that are result from matrix-vector 248 * products. 249 */ 250 std::vector<IndexSet> 251 locally_owned_range_indices() const; 252 253 /** 254 * Return a reference to the MPI communicator object in use with this 255 * matrix. 256 */ 257 const MPI_Comm & 258 get_mpi_communicator() const; 259 260 /** 261 * Make the clear() function in the base class visible, though it is 262 * protected. 263 */ 264 using BlockMatrixBase<SparseMatrix>::clear; 265 }; 266 267 268 269 /*@}*/ 270 271 // ------------- inline and template functions ----------------- 272 273 inline BlockSparseMatrix & 274 BlockSparseMatrix::operator=(const double d) 275 { 276 Assert(d == 0, ExcScalarAssignmentOnlyForZeroValue()); 277 278 for (size_type r = 0; r < this->n_block_rows(); ++r) 279 for (size_type c = 0; c < this->n_block_cols(); ++c) 280 this->block(r, c) = d; 281 282 return *this; 283 } 284 285 286 287 inline void vmult(BlockVector & dst,const BlockVector & src)288 BlockSparseMatrix::vmult(BlockVector &dst, const BlockVector &src) const 289 { 290 BaseClass::vmult_block_block(dst, src); 291 } 292 293 294 295 inline void vmult(BlockVector & dst,const Vector & src)296 BlockSparseMatrix::vmult(BlockVector &dst, const Vector &src) const 297 { 298 BaseClass::vmult_block_nonblock(dst, src); 299 } 300 301 302 303 inline void vmult(Vector & dst,const BlockVector & src)304 BlockSparseMatrix::vmult(Vector &dst, const BlockVector &src) const 305 { 306 BaseClass::vmult_nonblock_block(dst, src); 307 } 308 309 310 311 inline void vmult(Vector & dst,const Vector & src)312 BlockSparseMatrix::vmult(Vector &dst, const Vector &src) const 313 { 314 BaseClass::vmult_nonblock_nonblock(dst, src); 315 } 316 317 318 inline void Tvmult(BlockVector & dst,const BlockVector & src)319 BlockSparseMatrix::Tvmult(BlockVector &dst, const BlockVector &src) const 320 { 321 BaseClass::Tvmult_block_block(dst, src); 322 } 323 324 325 326 inline void Tvmult(BlockVector & dst,const Vector & src)327 BlockSparseMatrix::Tvmult(BlockVector &dst, const Vector &src) const 328 { 329 BaseClass::Tvmult_block_nonblock(dst, src); 330 } 331 332 333 334 inline void Tvmult(Vector & dst,const BlockVector & src)335 BlockSparseMatrix::Tvmult(Vector &dst, const BlockVector &src) const 336 { 337 BaseClass::Tvmult_nonblock_block(dst, src); 338 } 339 340 341 342 inline void Tvmult(Vector & dst,const Vector & src)343 BlockSparseMatrix::Tvmult(Vector &dst, const Vector &src) const 344 { 345 BaseClass::Tvmult_nonblock_nonblock(dst, src); 346 } 347 348 } // namespace MPI 349 350 } // namespace PETScWrappers 351 352 353 DEAL_II_NAMESPACE_CLOSE 354 355 356 #endif // DEAL_II_WITH_PETSC 357 358 #endif // dealii_petsc_block_sparse_matrix_h 359