1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2008 - 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_trilinos_block_sparse_matrix_h 17 #define dealii_trilinos_block_sparse_matrix_h 18 19 20 #include <deal.II/base/config.h> 21 22 #ifdef DEAL_II_WITH_TRILINOS 23 24 # include <deal.II/base/template_constraints.h> 25 26 # include <deal.II/lac/block_matrix_base.h> 27 # include <deal.II/lac/exceptions.h> 28 # include <deal.II/lac/full_matrix.h> 29 # include <deal.II/lac/trilinos_parallel_block_vector.h> 30 # include <deal.II/lac/trilinos_sparse_matrix.h> 31 32 # include <cmath> 33 34 DEAL_II_NAMESPACE_OPEN 35 36 // forward declarations 37 # ifndef DOXYGEN 38 class BlockSparsityPattern; 39 template <typename number> 40 class BlockSparseMatrix; 41 # endif 42 43 namespace TrilinosWrappers 44 { 45 /*! @addtogroup TrilinosWrappers 46 *@{ 47 */ 48 49 /** 50 * Blocked sparse matrix based on the TrilinosWrappers::SparseMatrix class. 51 * This class implements the functions that are specific to the Trilinos 52 * SparseMatrix base objects for a blocked sparse matrix, and leaves the 53 * actual work relaying most of the calls to the individual blocks to the 54 * functions implemented in the base class. See there also for a description 55 * of when this class is useful. 56 * 57 * In contrast to the deal.II-type SparseMatrix class, the Trilinos matrices 58 * do not have external objects for the sparsity patterns. Thus, one does 59 * not determine the size of the individual blocks of a block matrix of this 60 * type by attaching a block sparsity pattern, but by calling reinit() to 61 * set the number of blocks and then by setting the size of each block 62 * separately. In order to fix the data structures of the block matrix, it 63 * is then necessary to let it know that we have changed the sizes of the 64 * underlying matrices. For this, one has to call the collect_sizes() 65 * function, for much the same reason as is documented with the 66 * BlockSparsityPattern class. 67 * 68 * @ingroup Matrix1 @see 69 * @ref GlossBlockLA "Block (linear algebra)" 70 */ 71 class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix> 72 { 73 public: 74 /** 75 * Typedef the base class for simpler access to its own alias. 76 */ 77 using BaseClass = BlockMatrixBase<SparseMatrix>; 78 79 /** 80 * Typedef the type of the underlying matrix. 81 */ 82 using BlockType = BaseClass::BlockType; 83 84 /** 85 * Import the alias from the base class. 86 */ 87 using value_type = BaseClass::value_type; 88 using pointer = BaseClass::pointer; 89 using const_pointer = BaseClass::const_pointer; 90 using reference = BaseClass::reference; 91 using const_reference = BaseClass::const_reference; 92 using size_type = BaseClass::size_type; 93 using iterator = BaseClass::iterator; 94 using const_iterator = BaseClass::const_iterator; 95 96 /** 97 * Constructor; initializes the matrix to be empty, without any structure, 98 * i.e. the matrix is not usable at all. This constructor is therefore 99 * only useful for matrices which are members of a class. All other 100 * matrices should be created at a point in the data flow where all 101 * necessary information is available. 102 * 103 * You have to initialize the matrix before usage with 104 * reinit(BlockSparsityPattern). The number of blocks per row and column 105 * are then determined by that function. 106 */ 107 BlockSparseMatrix() = default; 108 109 /** 110 * Destructor. 111 */ 112 ~BlockSparseMatrix() override; 113 114 /** 115 * Pseudo copy operator only copying empty objects. The sizes of the block 116 * matrices need to be the same. 117 */ 118 BlockSparseMatrix & 119 operator=(const BlockSparseMatrix &) = default; 120 121 /** 122 * This operator assigns a scalar to a matrix. Since this does usually not 123 * make much sense (should we set all matrix entries to this value? Only 124 * the nonzero entries of the sparsity pattern?), this operation is only 125 * allowed if the actual value to be assigned is zero. This operator only 126 * exists to allow for the obvious notation <tt>matrix=0</tt>, which sets 127 * all elements of the matrix to zero, but keep the sparsity pattern 128 * previously used. 129 */ 130 BlockSparseMatrix & 131 operator=(const double d); 132 133 /** 134 * Resize the matrix, by setting the number of block rows and columns. 135 * This deletes all blocks and replaces them with uninitialized ones, i.e. 136 * ones for which also the sizes are not yet set. You have to do that by 137 * calling the @p reinit functions of the blocks themselves. Do not forget 138 * to call collect_sizes() after that on this object. 139 * 140 * The reason that you have to set sizes of the blocks yourself is that 141 * the sizes may be varying, the maximum number of elements per row may be 142 * varying, etc. It is simpler not to reproduce the interface of the @p 143 * SparsityPattern class here but rather let the user call whatever 144 * function they desire. 145 */ 146 void 147 reinit(const size_type n_block_rows, const size_type n_block_columns); 148 149 /** 150 * Resize the matrix, by using an array of index sets to determine the 151 * %parallel distribution of the individual matrices. This function 152 * assumes that a quadratic block matrix is generated. 153 */ 154 template <typename BlockSparsityPatternType> 155 void 156 reinit(const std::vector<IndexSet> & input_maps, 157 const BlockSparsityPatternType &block_sparsity_pattern, 158 const MPI_Comm & communicator = MPI_COMM_WORLD, 159 const bool exchange_data = false); 160 161 /** 162 * Resize the matrix and initialize it by the given sparsity pattern. 163 * Since no distribution map is given, the result is a block matrix for 164 * which all elements are stored locally. 165 */ 166 template <typename BlockSparsityPatternType> 167 void 168 reinit(const BlockSparsityPatternType &block_sparsity_pattern); 169 170 /** 171 * This function initializes the Trilinos matrix using the deal.II sparse 172 * matrix and the entries stored therein. It uses a threshold to copy only 173 * elements whose modulus is larger than the threshold (so zeros in the 174 * deal.II matrix can be filtered away). 175 */ 176 void 177 reinit( 178 const std::vector<IndexSet> & parallel_partitioning, 179 const ::dealii::BlockSparseMatrix<double> &dealii_block_sparse_matrix, 180 const MPI_Comm & communicator = MPI_COMM_WORLD, 181 const double drop_tolerance = 1e-13); 182 183 /** 184 * This function initializes the Trilinos matrix using the deal.II sparse 185 * matrix and the entries stored therein. It uses a threshold to copy only 186 * elements whose modulus is larger than the threshold (so zeros in the 187 * deal.II matrix can be filtered away). Since no Epetra_Map is given, all 188 * the elements will be locally stored. 189 */ 190 void 191 reinit(const ::dealii::BlockSparseMatrix<double> &deal_ii_sparse_matrix, 192 const double drop_tolerance = 1e-13); 193 194 /** 195 * Return the state of the matrix, i.e., whether compress() needs to be 196 * called after an operation requiring data exchange. Does only return 197 * non-true values when used in <tt>debug</tt> mode, since it is quite 198 * expensive to keep track of all operations that lead to the need for 199 * compress(). 200 */ 201 bool 202 is_compressed() const; 203 204 /** 205 * This function collects the sizes of the sub-objects and stores them in 206 * internal arrays, in order to be able to relay global indices into the 207 * matrix to indices into the subobjects. You *must* call this function 208 * each time after you have changed the size of the sub-objects. Note that 209 * this is a collective operation, i.e., it needs to be called on all MPI 210 * processes. This command internally calls the method 211 * <tt>compress()</tt>, so you don't need to call that function in case 212 * you use <tt>collect_sizes()</tt>. 213 */ 214 void 215 collect_sizes(); 216 217 /** 218 * Return the total number of nonzero elements of this matrix (summed 219 * over all MPI processes). 220 */ 221 size_type 222 n_nonzero_elements() const; 223 224 /** 225 * Return the MPI communicator object in use with this matrix. 226 */ 227 MPI_Comm 228 get_mpi_communicator() const; 229 230 /** 231 * Return the partitioning of the domain space for the individual blocks of 232 * this matrix, i.e., the partitioning of the block vectors this matrix has 233 * to be multiplied with. 234 */ 235 std::vector<IndexSet> 236 locally_owned_domain_indices() const; 237 238 /** 239 * Return the partitioning of the range space for the individual blocks of 240 * this matrix, i.e., the partitioning of the block vectors that result 241 * from matrix-vector products. 242 */ 243 std::vector<IndexSet> 244 locally_owned_range_indices() const; 245 246 /** 247 * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this 248 * matrix. The vector types can be block vectors or non-block vectors 249 * (only if the matrix has only one row or column, respectively), and need 250 * to define TrilinosWrappers::SparseMatrix::vmult. 251 */ 252 template <typename VectorType1, typename VectorType2> 253 void 254 vmult(VectorType1 &dst, const VectorType2 &src) const; 255 256 /** 257 * Matrix-vector multiplication: let $dst = M^T*src$ with $M$ being this 258 * matrix. This function does the same as vmult() but takes the transposed 259 * matrix. 260 */ 261 template <typename VectorType1, typename VectorType2> 262 void 263 Tvmult(VectorType1 &dst, const VectorType2 &src) const; 264 265 /** 266 * Compute the residual of an equation <i>Mx=b</i>, where the residual is 267 * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The 268 * <i>l<sub>2</sub></i> norm of the residual vector is returned. 269 * 270 * Source <i>x</i> and destination <i>dst</i> must not be the same vector. 271 * 272 * Note that both vectors have to be distributed vectors generated using 273 * the same Map as was used for the matrix. 274 * 275 * This function only applicable if the matrix only has one block row. 276 */ 277 TrilinosScalar 278 residual(MPI::BlockVector & dst, 279 const MPI::BlockVector &x, 280 const MPI::BlockVector &b) const; 281 282 /** 283 * Compute the residual of an equation <i>Mx=b</i>, where the residual is 284 * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The 285 * <i>l<sub>2</sub></i> norm of the residual vector is returned. 286 * 287 * This function is only applicable if the matrix only has one block row. 288 */ 289 TrilinosScalar 290 residual(MPI::BlockVector & dst, 291 const MPI::Vector & x, 292 const MPI::BlockVector &b) const; 293 294 /** 295 * Compute the residual of an equation <i>Mx=b</i>, where the residual is 296 * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The 297 * <i>l<sub>2</sub></i> norm of the residual vector is returned. 298 * 299 * This function is only applicable if the matrix only has one block column. 300 */ 301 TrilinosScalar 302 residual(MPI::Vector & dst, 303 const MPI::BlockVector &x, 304 const MPI::Vector & b) const; 305 306 /** 307 * Compute the residual of an equation <i>Mx=b</i>, where the residual is 308 * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The 309 * <i>l<sub>2</sub></i> norm of the residual vector is returned. 310 * 311 * This function is only applicable if the matrix only has one block. 312 */ 313 TrilinosScalar 314 residual(MPI::Vector & dst, 315 const MPI::Vector &x, 316 const MPI::Vector &b) const; 317 318 /** 319 * Make the clear() function in the base class visible, though it is 320 * protected. 321 */ 322 using BlockMatrixBase<SparseMatrix>::clear; 323 324 /** 325 * @addtogroup Exceptions 326 * @{ 327 */ 328 329 /** 330 * Exception 331 */ 332 DeclException4(ExcIncompatibleRowNumbers, 333 int, 334 int, 335 int, 336 int, 337 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3 338 << ',' << arg4 << "] have differing row numbers."); 339 340 /** 341 * Exception 342 */ 343 DeclException4(ExcIncompatibleColNumbers, 344 int, 345 int, 346 int, 347 int, 348 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3 349 << ',' << arg4 << "] have differing column numbers."); 350 ///@} 351 352 private: 353 /** 354 * Internal version of (T)vmult with two block vectors 355 */ 356 template <typename VectorType1, typename VectorType2> 357 void 358 vmult(VectorType1 & dst, 359 const VectorType2 &src, 360 const bool transpose, 361 const std::integral_constant<bool, true>, 362 const std::integral_constant<bool, true>) const; 363 364 /** 365 * Internal version of (T)vmult where the source vector is a block vector 366 * but the destination vector is a non-block vector 367 */ 368 template <typename VectorType1, typename VectorType2> 369 void 370 vmult(VectorType1 & dst, 371 const VectorType2 &src, 372 const bool transpose, 373 const std::integral_constant<bool, false>, 374 const std::integral_constant<bool, true>) const; 375 376 /** 377 * Internal version of (T)vmult where the source vector is a non-block 378 * vector but the destination vector is a block vector 379 */ 380 template <typename VectorType1, typename VectorType2> 381 void 382 vmult(VectorType1 & dst, 383 const VectorType2 &src, 384 const bool transpose, 385 const std::integral_constant<bool, true>, 386 const std::integral_constant<bool, false>) const; 387 388 /** 389 * Internal version of (T)vmult where both source vector and the 390 * destination vector are non-block vectors (only defined if the matrix 391 * consists of only one block) 392 */ 393 template <typename VectorType1, typename VectorType2> 394 void 395 vmult(VectorType1 & dst, 396 const VectorType2 &src, 397 const bool transpose, 398 const std::integral_constant<bool, false>, 399 const std::integral_constant<bool, false>) const; 400 }; 401 402 403 404 /*@}*/ 405 406 // ------------- inline and template functions ----------------- 407 408 409 410 inline BlockSparseMatrix & 411 BlockSparseMatrix::operator=(const double d) 412 { 413 Assert(d == 0, ExcScalarAssignmentOnlyForZeroValue()); 414 415 for (size_type r = 0; r < this->n_block_rows(); ++r) 416 for (size_type c = 0; c < this->n_block_cols(); ++c) 417 this->block(r, c) = d; 418 419 return *this; 420 } 421 422 423 424 inline bool is_compressed()425 BlockSparseMatrix::is_compressed() const 426 { 427 bool compressed = true; 428 for (size_type row = 0; row < n_block_rows(); ++row) 429 for (size_type col = 0; col < n_block_cols(); ++col) 430 if (block(row, col).is_compressed() == false) 431 { 432 compressed = false; 433 break; 434 } 435 436 return compressed; 437 } 438 439 440 441 template <typename VectorType1, typename VectorType2> 442 inline void vmult(VectorType1 & dst,const VectorType2 & src)443 BlockSparseMatrix::vmult(VectorType1 &dst, const VectorType2 &src) const 444 { 445 vmult(dst, 446 src, 447 false, 448 std::integral_constant<bool, IsBlockVector<VectorType1>::value>(), 449 std::integral_constant<bool, IsBlockVector<VectorType2>::value>()); 450 } 451 452 453 454 template <typename VectorType1, typename VectorType2> 455 inline void Tvmult(VectorType1 & dst,const VectorType2 & src)456 BlockSparseMatrix::Tvmult(VectorType1 &dst, const VectorType2 &src) const 457 { 458 vmult(dst, 459 src, 460 true, 461 std::integral_constant<bool, IsBlockVector<VectorType1>::value>(), 462 std::integral_constant<bool, IsBlockVector<VectorType2>::value>()); 463 } 464 465 466 467 template <typename VectorType1, typename VectorType2> 468 inline void vmult(VectorType1 & dst,const VectorType2 & src,const bool transpose,std::integral_constant<bool,true>,std::integral_constant<bool,true>)469 BlockSparseMatrix::vmult(VectorType1 & dst, 470 const VectorType2 &src, 471 const bool transpose, 472 std::integral_constant<bool, true>, 473 std::integral_constant<bool, true>) const 474 { 475 if (transpose == true) 476 BaseClass::Tvmult_block_block(dst, src); 477 else 478 BaseClass::vmult_block_block(dst, src); 479 } 480 481 482 483 template <typename VectorType1, typename VectorType2> 484 inline void vmult(VectorType1 & dst,const VectorType2 & src,const bool transpose,std::integral_constant<bool,false>,std::integral_constant<bool,true>)485 BlockSparseMatrix::vmult(VectorType1 & dst, 486 const VectorType2 &src, 487 const bool transpose, 488 std::integral_constant<bool, false>, 489 std::integral_constant<bool, true>) const 490 { 491 if (transpose == true) 492 BaseClass::Tvmult_nonblock_block(dst, src); 493 else 494 BaseClass::vmult_nonblock_block(dst, src); 495 } 496 497 498 499 template <typename VectorType1, typename VectorType2> 500 inline void vmult(VectorType1 & dst,const VectorType2 & src,const bool transpose,std::integral_constant<bool,true>,std::integral_constant<bool,false>)501 BlockSparseMatrix::vmult(VectorType1 & dst, 502 const VectorType2 &src, 503 const bool transpose, 504 std::integral_constant<bool, true>, 505 std::integral_constant<bool, false>) const 506 { 507 if (transpose == true) 508 BaseClass::Tvmult_block_nonblock(dst, src); 509 else 510 BaseClass::vmult_block_nonblock(dst, src); 511 } 512 513 514 515 template <typename VectorType1, typename VectorType2> 516 inline void vmult(VectorType1 & dst,const VectorType2 & src,const bool transpose,std::integral_constant<bool,false>,std::integral_constant<bool,false>)517 BlockSparseMatrix::vmult(VectorType1 & dst, 518 const VectorType2 &src, 519 const bool transpose, 520 std::integral_constant<bool, false>, 521 std::integral_constant<bool, false>) const 522 { 523 if (transpose == true) 524 BaseClass::Tvmult_nonblock_nonblock(dst, src); 525 else 526 BaseClass::vmult_nonblock_nonblock(dst, src); 527 } 528 529 530 531 inline std::vector<IndexSet> locally_owned_domain_indices()532 BlockSparseMatrix::locally_owned_domain_indices() const 533 { 534 Assert(this->n_block_cols() != 0, ExcNotInitialized()); 535 Assert(this->n_block_rows() != 0, ExcNotInitialized()); 536 537 std::vector<IndexSet> domain_indices; 538 for (size_type c = 0; c < this->n_block_cols(); ++c) 539 domain_indices.push_back( 540 this->sub_objects[0][c]->locally_owned_domain_indices()); 541 542 return domain_indices; 543 } 544 545 546 547 inline std::vector<IndexSet> locally_owned_range_indices()548 BlockSparseMatrix::locally_owned_range_indices() const 549 { 550 Assert(this->n_block_cols() != 0, ExcNotInitialized()); 551 Assert(this->n_block_rows() != 0, ExcNotInitialized()); 552 553 std::vector<IndexSet> range_indices; 554 for (size_type r = 0; r < this->n_block_rows(); ++r) 555 range_indices.push_back( 556 this->sub_objects[r][0]->locally_owned_range_indices()); 557 558 return range_indices; 559 } 560 561 562 563 namespace internal 564 { 565 namespace BlockLinearOperatorImplementation 566 { 567 /** 568 * This is an extension class to BlockLinearOperators for Trilinos block 569 * sparse matrices. 570 * 571 * @note This class does very little at the moment other than to check 572 * that the correct Payload type for each subblock has been chosen 573 * correctly. Further extensions to the class may be necessary in the 574 * future in order to add further functionality to BlockLinearOperators 575 * while retaining compatibility with the Trilinos sparse matrix and 576 * preconditioner classes. 577 * 578 * 579 * @ingroup TrilinosWrappers 580 */ 581 template <typename PayloadBlockType> 582 class TrilinosBlockPayload 583 { 584 public: 585 /** 586 * Type of payload held by each subblock 587 */ 588 using BlockType = PayloadBlockType; 589 590 /** 591 * Default constructor 592 * 593 * This simply checks that the payload for each block has been chosen 594 * correctly (i.e. is of type TrilinosPayload). Apart from this, this 595 * class does not do anything in particular and needs no special 596 * configuration, we have only one generic constructor that can be 597 * called under any conditions. 598 */ 599 template <typename... Args> TrilinosBlockPayload(const Args &...)600 TrilinosBlockPayload(const Args &...) 601 { 602 static_assert( 603 std::is_same< 604 PayloadBlockType, 605 internal::LinearOperatorImplementation::TrilinosPayload>::value, 606 "TrilinosBlockPayload can only accept a payload of type TrilinosPayload."); 607 } 608 }; 609 610 } // namespace BlockLinearOperatorImplementation 611 } /* namespace internal */ 612 613 614 } /* namespace TrilinosWrappers */ 615 616 617 DEAL_II_NAMESPACE_CLOSE 618 619 #endif // DEAL_II_WITH_TRILINOS 620 621 #endif // dealii_trilinos_block_sparse_matrix_h 622