00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__block_sparse_matrix_h
00013 #define __deal2__block_sparse_matrix_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/table.h>
00018 #include <deal.II/lac/block_matrix_base.h>
00019 #include <deal.II/lac/block_vector.h>
00020 #include <deal.II/lac/sparse_matrix.h>
00021 #include <deal.II/lac/block_sparsity_pattern.h>
00022 #include <deal.II/lac/exceptions.h>
00023
00024 #include <cmath>
00025
00026 DEAL_II_NAMESPACE_OPEN
00027
00028
00044 template <typename number>
00045 class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix<number> >
00046 {
00047 public:
00052 typedef BlockMatrixBase<SparseMatrix<number> > BaseClass;
00053
00058 typedef typename BaseClass::BlockType BlockType;
00059
00064 typedef typename BaseClass::value_type value_type;
00065 typedef typename BaseClass::pointer pointer;
00066 typedef typename BaseClass::const_pointer const_pointer;
00067 typedef typename BaseClass::reference reference;
00068 typedef typename BaseClass::const_reference const_reference;
00069 typedef typename BaseClass::size_type size_type;
00070 typedef typename BaseClass::iterator iterator;
00071 typedef typename BaseClass::const_iterator const_iterator;
00072
00098 BlockSparseMatrix ();
00099
00121 BlockSparseMatrix (const BlockSparsityPattern &sparsity);
00122
00126 virtual ~BlockSparseMatrix ();
00127
00128
00129
00135 BlockSparseMatrix &
00136 operator = (const BlockSparseMatrix &);
00137
00153 BlockSparseMatrix &
00154 operator = (const double d);
00155
00168 void clear ();
00169
00187 virtual void reinit (const BlockSparsityPattern &sparsity);
00189
00201 bool empty () const;
00202
00207 unsigned int get_row_length (const unsigned int row) const;
00208
00218 unsigned int n_nonzero_elements () const;
00219
00227 unsigned int n_actually_nonzero_elements (const double threshold = 0.0) const;
00228
00241 const BlockSparsityPattern &
00242 get_sparsity_pattern () const;
00243
00249 std::size_t memory_consumption () const;
00251
00261 template <typename block_number>
00262 void vmult (BlockVector<block_number> &dst,
00263 const BlockVector<block_number> &src) const;
00264
00272 template <typename block_number,
00273 typename nonblock_number>
00274 void vmult (BlockVector<block_number> &dst,
00275 const Vector<nonblock_number> &src) const;
00276
00284 template <typename block_number,
00285 typename nonblock_number>
00286 void vmult (Vector<nonblock_number> &dst,
00287 const BlockVector<block_number> &src) const;
00288
00296 template <typename nonblock_number>
00297 void vmult (Vector<nonblock_number> &dst,
00298 const Vector<nonblock_number> &src) const;
00299
00308 template <typename block_number>
00309 void Tvmult (BlockVector<block_number> &dst,
00310 const BlockVector<block_number> &src) const;
00311
00319 template <typename block_number,
00320 typename nonblock_number>
00321 void Tvmult (BlockVector<block_number> &dst,
00322 const Vector<nonblock_number> &src) const;
00323
00331 template <typename block_number,
00332 typename nonblock_number>
00333 void Tvmult (Vector<nonblock_number> &dst,
00334 const BlockVector<block_number> &src) const;
00335
00343 template <typename nonblock_number>
00344 void Tvmult (Vector<nonblock_number> &dst,
00345 const Vector<nonblock_number> &src) const;
00347
00367 template <class BlockVectorType>
00368 void precondition_Jacobi (BlockVectorType &dst,
00369 const BlockVectorType &src,
00370 const number omega = 1.) const;
00371
00379 template <typename number2>
00380 void precondition_Jacobi (Vector<number2> &dst,
00381 const Vector<number2> &src,
00382 const number omega = 1.) const;
00384
00430 void print_formatted (std::ostream &out,
00431 const unsigned int precision = 3,
00432 const bool scientific = true,
00433 const unsigned int width = 0,
00434 const char *zero_string = " ",
00435 const double denominator = 1.) const;
00437
00443 DeclException0 (ExcBlockDimensionMismatch);
00445
00446 private:
00456 SmartPointer<const BlockSparsityPattern,BlockSparseMatrix<number> > sparsity_pattern;
00457 };
00458
00459
00460
00462
00463
00464
00465
00466 template <typename number>
00467 inline
00468 BlockSparseMatrix<number> &
00469 BlockSparseMatrix<number>::operator = (const double d)
00470 {
00471 Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
00472
00473 for (unsigned int r=0; r<this->n_block_rows(); ++r)
00474 for (unsigned int c=0; c<this->n_block_cols(); ++c)
00475 this->block(r,c) = d;
00476
00477 return *this;
00478 }
00479
00480
00481
00482 template <typename number>
00483 template <typename block_number>
00484 inline
00485 void
00486 BlockSparseMatrix<number>::vmult (BlockVector<block_number> &dst,
00487 const BlockVector<block_number> &src) const
00488 {
00489 BaseClass::vmult_block_block (dst, src);
00490 }
00491
00492
00493
00494 template <typename number>
00495 template <typename block_number,
00496 typename nonblock_number>
00497 inline
00498 void
00499 BlockSparseMatrix<number>::vmult (BlockVector<block_number> &dst,
00500 const Vector<nonblock_number> &src) const
00501 {
00502 BaseClass::vmult_block_nonblock (dst, src);
00503 }
00504
00505
00506
00507 template <typename number>
00508 template <typename block_number,
00509 typename nonblock_number>
00510 inline
00511 void
00512 BlockSparseMatrix<number>::vmult (Vector<nonblock_number> &dst,
00513 const BlockVector<block_number> &src) const
00514 {
00515 BaseClass::vmult_nonblock_block (dst, src);
00516 }
00517
00518
00519
00520 template <typename number>
00521 template <typename nonblock_number>
00522 inline
00523 void
00524 BlockSparseMatrix<number>::vmult (Vector<nonblock_number> &dst,
00525 const Vector<nonblock_number> &src) const
00526 {
00527 BaseClass::vmult_nonblock_nonblock (dst, src);
00528 }
00529
00530
00531
00532 template <typename number>
00533 template <typename block_number>
00534 inline
00535 void
00536 BlockSparseMatrix<number>::Tvmult (BlockVector<block_number> &dst,
00537 const BlockVector<block_number> &src) const
00538 {
00539 BaseClass::Tvmult_block_block (dst, src);
00540 }
00541
00542
00543
00544 template <typename number>
00545 template <typename block_number,
00546 typename nonblock_number>
00547 inline
00548 void
00549 BlockSparseMatrix<number>::Tvmult (BlockVector<block_number> &dst,
00550 const Vector<nonblock_number> &src) const
00551 {
00552 BaseClass::Tvmult_block_nonblock (dst, src);
00553 }
00554
00555
00556
00557 template <typename number>
00558 template <typename block_number,
00559 typename nonblock_number>
00560 inline
00561 void
00562 BlockSparseMatrix<number>::Tvmult (Vector<nonblock_number> &dst,
00563 const BlockVector<block_number> &src) const
00564 {
00565 BaseClass::Tvmult_nonblock_block (dst, src);
00566 }
00567
00568
00569
00570 template <typename number>
00571 template <typename nonblock_number>
00572 inline
00573 void
00574 BlockSparseMatrix<number>::Tvmult (Vector<nonblock_number> &dst,
00575 const Vector<nonblock_number> &src) const
00576 {
00577 BaseClass::Tvmult_nonblock_nonblock (dst, src);
00578 }
00579
00580
00581
00582 template <typename number>
00583 template <class BlockVectorType>
00584 inline
00585 void
00586 BlockSparseMatrix<number>::
00587 precondition_Jacobi (BlockVectorType &dst,
00588 const BlockVectorType &src,
00589 const number omega) const
00590 {
00591 Assert (this->n_block_rows() == this->n_block_cols(), ExcNotQuadratic());
00592 Assert (dst.n_blocks() == this->n_block_rows(),
00593 ExcDimensionMismatch(dst.n_blocks(), this->n_block_rows()));
00594 Assert (src.n_blocks() == this->n_block_cols(),
00595 ExcDimensionMismatch(src.n_blocks(), this->n_block_cols()));
00596
00597
00598
00599 for (unsigned int i=0; i<this->n_block_rows(); ++i)
00600 this->block(i,i).precondition_Jacobi (dst.block(i),
00601 src.block(i),
00602 omega);
00603 }
00604
00605
00606
00607 template <typename number>
00608 template <typename number2>
00609 inline
00610 void
00611 BlockSparseMatrix<number>::
00612 precondition_Jacobi (Vector<number2> &dst,
00613 const Vector<number2> &src,
00614 const number omega) const
00615 {
00616
00617
00618
00619 Assert (this->n_block_cols() == 1,
00620 ExcMessage ("This function only works if the matrix has "
00621 "a single block"));
00622 Assert (this->n_block_rows() == 1,
00623 ExcMessage ("This function only works if the matrix has "
00624 "a single block"));
00625
00626
00627
00628 this->block(0,0).precondition_Jacobi (dst, src, omega);
00629 }
00630
00631
00632 DEAL_II_NAMESPACE_CLOSE
00633
00634 #endif // __deal2__block_sparse_matrix_h
00635