00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__block_matrix_array_h
00013 #define __deal2__block_matrix_array_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/subscriptor.h>
00017 #include <deal.II/base/table.h>
00018
00019 #include <deal.II/lac/pointer_matrix.h>
00020 #include <deal.II/lac/vector_memory.h>
00021
00022 #include <vector>
00023 #include <map>
00024 #include <string>
00025 #include <memory>
00026 #include <sstream>
00027
00028 DEAL_II_NAMESPACE_OPEN
00029
00030 template <typename> class BlockVector;
00031 template <typename> class Vector;
00032
00122 template <typename number = double>
00123 class BlockMatrixArray : public Subscriptor
00124 {
00125 public:
00132 BlockMatrixArray ();
00133
00138 BlockMatrixArray (const unsigned int n_block_rows,
00139 const unsigned int n_block_cols);
00140
00151 void initialize (const unsigned int n_block_rows,
00152 const unsigned int n_block_cols);
00153
00161 BlockMatrixArray (const unsigned int n_block_rows,
00162 const unsigned int n_block_cols,
00163 VectorMemory<Vector<number> >& mem);
00164
00175 void initialize (const unsigned int n_block_rows,
00176 const unsigned int n_block_cols,
00177 VectorMemory<Vector<number> >& mem);
00178
00183 void reinit (const unsigned int n_block_rows,
00184 const unsigned int n_block_cols);
00185
00205 template <class MATRIX>
00206 void enter (const MATRIX &matrix,
00207 const unsigned int row,
00208 const unsigned int col,
00209 const double prefix = 1.,
00210 const bool transpose = false);
00211
00221 template <class MATRIX>
00222 void enter_aux (VectorMemory<Vector<number> >& mem,
00223 const MATRIX &matrix,
00224 const unsigned int row,
00225 const unsigned int col,
00226 const double prefix = 1.,
00227 const bool transpose = false);
00228
00229
00234 void clear();
00235
00240 unsigned int n_block_rows () const;
00241
00246 unsigned int n_block_cols () const;
00247
00251 void vmult (BlockVector<number>& dst,
00252 const BlockVector<number>& src) const;
00253
00258 void vmult_add (BlockVector<number>& dst,
00259 const BlockVector<number>& src) const;
00260
00265 void Tvmult (BlockVector<number>& dst,
00266 const BlockVector<number>& src) const;
00267
00273 void Tvmult_add (BlockVector<number>& dst,
00274 const BlockVector<number>& src) const;
00275
00281 number matrix_scalar_product (const BlockVector<number>& u,
00282 const BlockVector<number>& v) const;
00283
00290 number matrix_norm_square (const BlockVector<number>& u) const;
00291
00309 template <class STREAM>
00310 void print_latex (STREAM& out) const;
00311
00312 protected:
00325 class Entry
00326 {
00327 public:
00335 template<class MATRIX>
00336 Entry (const MATRIX& matrix,
00337 unsigned row, unsigned int col,
00338 double prefix, bool transpose);
00339
00353 Entry(const Entry&);
00354
00361 ~Entry();
00362
00367 unsigned int row;
00368
00373 unsigned int col;
00374
00379 double prefix;
00380
00386 bool transpose;
00387
00391 PointerMatrixBase<Vector<number> >* matrix;
00392 };
00393
00398 std::vector<Entry> entries;
00399
00400 private:
00404 unsigned int block_rows;
00408 unsigned int block_cols;
00409 };
00410
00471 template <typename number = double>
00472 class BlockTrianglePrecondition
00473 : private BlockMatrixArray<number>
00474 {
00475 public:
00482 BlockTrianglePrecondition ();
00483
00490 BlockTrianglePrecondition (unsigned int n_blocks);
00491
00501 BlockTrianglePrecondition (unsigned int n_block_rows,
00502 VectorMemory<Vector<number> >& mem,
00503 bool backward = false);
00504
00515 void initialize (const unsigned int n_block_rows,
00516 VectorMemory<Vector<number> >& mem,
00517 bool backward = false);
00518
00523 void reinit(const unsigned int n_block_rows);
00524
00525
00533 template <class MATRIX>
00534 void enter (const MATRIX &matrix,
00535 const unsigned int row,
00536 const unsigned int col,
00537 const double prefix = 1.,
00538 const bool transpose = false);
00539
00551 template <class MATRIX>
00552 void enter_aux (VectorMemory<Vector<double> >& mem,
00553 const MATRIX &matrix,
00554 const unsigned int row,
00555 const unsigned int col,
00556 const double prefix = 1.,
00557 const bool transpose = false);
00558
00562 void vmult (BlockVector<number>& dst,
00563 const BlockVector<number>& src) const;
00564
00569 void vmult_add (BlockVector<number>& dst,
00570 const BlockVector<number>& src) const;
00571
00575 void Tvmult (BlockVector<number>& dst,
00576 const BlockVector<number>& src) const;
00577
00582 void Tvmult_add (BlockVector<number>& dst,
00583 const BlockVector<number>& src) const;
00584
00588 using BlockMatrixArray<number>::print_latex;
00589
00593 using BlockMatrixArray<number>::n_block_rows;
00594
00598 using BlockMatrixArray<number>::n_block_cols;
00599 using BlockMatrixArray<number>::clear;
00600 using BlockMatrixArray<number>::Subscriptor::subscribe;
00601 using BlockMatrixArray<number>::Subscriptor::unsubscribe;
00602
00613 DeclException1(ExcNoDiagonal,
00614 unsigned int,
00615 << "No diagonal entry was added for block " << arg1);
00616
00624 DeclException1(ExcMultipleDiagonal,
00625 unsigned int,
00626 << "Inverse diagonal entries may not be added in block "
00627 << arg1);
00629 private:
00636 void do_row (BlockVector<number>& dst,
00637 unsigned int row_num) const;
00638
00642 bool backward;
00643 };
00644
00645
00646 #ifndef DOXYGEN
00647
00648
00649 template <typename number>
00650 template <class MATRIX>
00651 inline
00652 BlockMatrixArray<number>::Entry::Entry (
00653 const MATRIX& m,
00654 unsigned int row,
00655 unsigned int col,
00656 double prefix,
00657 bool transpose)
00658 :
00659 row (row),
00660 col (col),
00661 prefix (prefix),
00662 transpose (transpose),
00663 matrix (new_pointer_matrix_base(m, Vector<number>(), typeid(*this).name()))
00664 {}
00665
00666
00667
00668 template <typename number>
00669 template <class MATRIX>
00670 inline
00671 void
00672 BlockMatrixArray<number>::enter (
00673 const MATRIX& matrix,
00674 unsigned int row,
00675 unsigned int col,
00676 double prefix,
00677 bool transpose)
00678 {
00679 Assert(row<n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
00680 Assert(col<n_block_cols(), ExcIndexRange(col, 0, n_block_cols()));
00681 entries.push_back(Entry(matrix, row, col, prefix, transpose));
00682 }
00683
00684
00685 template <typename number>
00686 template <class MATRIX>
00687 inline
00688 void
00689 BlockMatrixArray<number>::enter_aux (
00690 VectorMemory<Vector<number> >& mem,
00691 const MATRIX& matrix,
00692 unsigned int row,
00693 unsigned int col,
00694 double prefix,
00695 bool transpose)
00696 {
00697 Assert(row<n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
00698 Assert(col<n_block_cols(), ExcIndexRange(col, 0, n_block_cols()));
00699 entries.push_back(Entry(matrix, row, col, prefix, transpose, mem));
00700 }
00701
00702
00703 template<typename number>
00704 struct BlockMatrixArrayPointerMatrixLess
00705 {
00706 bool operator () (const PointerMatrixBase<Vector<number> >*a,
00707 const PointerMatrixBase<Vector<number> >* b) const
00708 {
00709 return *a < *b;
00710 }
00711 };
00712
00713
00714 template <typename number>
00715 template <class STREAM>
00716 inline
00717 void
00718 BlockMatrixArray<number>::print_latex (STREAM& out) const
00719 {
00720 out << "\\begin{array}{"
00721 << std::string(n_block_cols(), 'c')
00722 << "}" << std::endl;
00723
00724 Table<2,std::string> array(n_block_rows(), n_block_cols());
00725
00726 typedef std::map<const PointerMatrixBase<Vector<number> >*,
00727 std::string, BlockMatrixArrayPointerMatrixLess<number> > NameMap;
00728 NameMap matrix_names;
00729
00730 typename std::vector<Entry>::const_iterator m = entries.begin();
00731 typename std::vector<Entry>::const_iterator end = entries.end();
00732
00733 unsigned int matrix_number = 0;
00734 for (; m != end ; ++m)
00735 {
00736 if (matrix_names.find(m->matrix) == matrix_names.end())
00737 {
00738 std::pair<typename NameMap::iterator, bool> x =
00739 matrix_names.insert(
00740 std::pair<const PointerMatrixBase<Vector<number> >*, std::string> (m->matrix,
00741 std::string("M")));
00742 std::ostringstream stream;
00743 stream << matrix_number++;
00744
00745 x.first->second += stream.str();
00746 }
00747
00748 std::ostringstream stream;
00749
00750 if (array(m->row, m->col) != "" && m->prefix >= 0)
00751 stream << "+";
00752 if (m->prefix != 1.)
00753 stream << m->prefix << 'x';
00754 stream << matrix_names.find(m->matrix)->second;
00755
00756 if (m->transpose)
00757 stream << "^T";
00758
00759 array(m->row, m->col) += stream.str();
00760 }
00761 for (unsigned int i=0;i<n_block_rows();++i)
00762 for (unsigned int j=0;j<n_block_cols();++j)
00763 {
00764 out << '\t' << array(i,j);
00765 if (j==n_block_cols()-1)
00766 out << "\\\\" << std::endl;
00767 else
00768 out << " &";
00769 }
00770 out << "\\end{array}" << std::endl;
00771 }
00772
00773 template <typename number>
00774 template <class MATRIX>
00775 inline
00776 void
00777 BlockTrianglePrecondition<number>::enter (const MATRIX& matrix,
00778 unsigned row, unsigned int col,
00779 double prefix, bool transpose)
00780 {
00781 BlockMatrixArray<number>::enter(matrix, row, col, prefix, transpose);
00782 }
00783
00784
00785
00786 template <typename number>
00787 template <class MATRIX>
00788 inline
00789 void
00790 BlockTrianglePrecondition<number>::enter_aux (
00791 VectorMemory<Vector<double> >& mem,
00792 const MATRIX& matrix,
00793 unsigned int row,
00794 unsigned int col,
00795 double prefix,
00796 bool transpose)
00797 {
00798 BlockMatrixArray<number>::enter_aux(mem, matrix, row, col, prefix, transpose);
00799 }
00800
00801
00802
00803 #endif // DOXYGEN
00804
00805 DEAL_II_NAMESPACE_CLOSE
00806
00807 #endif
00808