00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__matrix_block_h
00013 #define __deal2__matrix_block_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/named_data.h>
00017 #include <deal.II/base/smartpointer.h>
00018 #include <deal.II/base/std_cxx1x/shared_ptr.h>
00019 #include <deal.II/base/memory_consumption.h>
00020 #include <deal.II/base/mg_level_object.h>
00021 #include <deal.II/lac/block_indices.h>
00022 #include <deal.II/lac/block_sparsity_pattern.h>
00023 #include <deal.II/lac/sparse_matrix.h>
00024 #include <deal.II/lac/full_matrix.h>
00025
00026
00027 DEAL_II_NAMESPACE_OPEN
00028
00029 template <class MATRIX> class MatrixBlock;
00030
00031 namespace internal
00032 {
00033 template <class MATRIX>
00034 void
00035 reinit(MatrixBlock<MATRIX>& v, const BlockSparsityPattern& p);
00036
00037 template <typename number>
00038 void
00039 reinit(MatrixBlock<::SparseMatrix<number> >&, const BlockSparsityPattern&);
00040 }
00041
00099 template <class MATRIX>
00100 class MatrixBlock
00101 : public Subscriptor
00102 {
00103 public:
00108 MatrixBlock();
00109
00113 MatrixBlock(const MatrixBlock<MATRIX>& M);
00114
00121 MatrixBlock(unsigned int i, unsigned int j);
00122
00134 void reinit(const BlockSparsityPattern& sparsity);
00135
00136 operator MATRIX&();
00137 operator const MATRIX&() const;
00138
00146 void add (const unsigned int i,
00147 const unsigned int j,
00148 const typename MATRIX::value_type value);
00149
00180 template <typename number>
00181 void add (const std::vector<unsigned int>& indices,
00182 const FullMatrix<number>& full_matrix,
00183 const bool elide_zero_values = true);
00184
00214 template <typename number>
00215 void add (const std::vector<unsigned int>& row_indices,
00216 const std::vector<unsigned int>& col_indices,
00217 const FullMatrix<number>& full_matrix,
00218 const bool elide_zero_values = true);
00219
00254 template <typename number>
00255 void add (const unsigned int row_index,
00256 const std::vector<unsigned int>& col_indices,
00257 const std::vector<number>& values,
00258 const bool elide_zero_values = true);
00259
00277 template <typename number>
00278 void add (const unsigned int row,
00279 const unsigned int n_cols,
00280 const unsigned int *col_indices,
00281 const number *values,
00282 const bool elide_zero_values = true,
00283 const bool col_indices_are_sorted = false);
00284
00293 template<class VECTOR>
00294 void vmult (VECTOR& w, const VECTOR& v) const;
00295
00304 template<class VECTOR>
00305 void vmult_add (VECTOR& w, const VECTOR& v) const;
00306
00315 template<class VECTOR>
00316 void Tvmult (VECTOR& w, const VECTOR& v) const;
00317
00326 template<class VECTOR>
00327 void Tvmult_add (VECTOR& w, const VECTOR& v) const;
00328
00332 std::size_t memory_consumption () const;
00333
00341 DeclException2(ExcBlockIndexMismatch, unsigned int, unsigned int,
00342 << "Block index " << arg1 << " does not match " << arg2);
00343
00349 unsigned int row;
00356 unsigned int column;
00357
00361 MATRIX matrix;
00362 private:
00370 BlockIndices row_indices;
00379 BlockIndices column_indices;
00380
00381 friend void internal::reinit<>(MatrixBlock<MATRIX>&, const BlockSparsityPattern&);
00382 };
00383
00384
00394 template <class MATRIX>
00395 class MatrixBlockVector
00396 :
00397 private NamedData<std_cxx1x::shared_ptr<MatrixBlock<MATRIX> > >
00398 {
00399 public:
00403 typedef MatrixBlock<MATRIX> value_type;
00404
00410 void add(unsigned int row, unsigned int column, const std::string& name);
00411
00419 void reinit(const BlockSparsityPattern& sparsity);
00420
00438 void clear (bool really_clean = false);
00439
00443 std::size_t memory_consumption () const;
00444
00449 const value_type& block(unsigned int i) const;
00450
00455 value_type& block(unsigned int i);
00456
00462 MATRIX& matrix(unsigned int i);
00463
00467 using NamedData<std_cxx1x::shared_ptr<value_type> >::subscribe;
00468 using NamedData<std_cxx1x::shared_ptr<value_type> >::unsubscribe;
00469 using NamedData<std_cxx1x::shared_ptr<value_type> >::size;
00470 using NamedData<std_cxx1x::shared_ptr<value_type> >::name;
00471 };
00472
00473
00483 template <class MATRIX>
00484 class MGMatrixBlockVector
00485 : public Subscriptor
00486 {
00487 public:
00491 typedef MGLevelObject<MatrixBlock<MATRIX> > value_type;
00507 MGMatrixBlockVector(const bool edge_matrices = false,
00508 const bool edge_flux_matrices = false);
00509
00513 unsigned int size () const;
00514
00523 void add(unsigned int row, unsigned int column,
00524 const std::string& name);
00525
00536 void reinit_matrix(const MGLevelObject<BlockSparsityPattern>& sparsity);
00548 void reinit_edge(const MGLevelObject<BlockSparsityPattern>& sparsity);
00560 void reinit_edge_flux(const MGLevelObject<BlockSparsityPattern>& sparsity);
00561
00579 void clear (bool really_clean = false);
00580
00585 const value_type& block(unsigned int i) const;
00586
00591 value_type& block(unsigned int i);
00592
00597 const value_type& block_in(unsigned int i) const;
00598
00603 value_type& block_in(unsigned int i);
00604
00609 const value_type& block_out(unsigned int i) const;
00610
00615 value_type& block_out(unsigned int i);
00616
00621 const value_type& block_up(unsigned int i) const;
00622
00627 value_type& block_up(unsigned int i);
00628
00633 const value_type& block_down(unsigned int i) const;
00634
00639 value_type& block_down(unsigned int i);
00640
00644 std::size_t memory_consumption () const;
00645 private:
00647 void clear_object(NamedData<MGLevelObject<MatrixBlock<MATRIX> > >&);
00648
00650 const bool edge_matrices;
00651
00653 const bool edge_flux_matrices;
00654
00656 NamedData<MGLevelObject<MatrixBlock<MATRIX> > > matrices;
00658 NamedData<MGLevelObject<MatrixBlock<MATRIX> > > matrices_in;
00660 NamedData<MGLevelObject<MatrixBlock<MATRIX> > > matrices_out;
00662 NamedData<MGLevelObject<MatrixBlock<MATRIX> > > flux_matrices_down;
00664 NamedData<MGLevelObject<MatrixBlock<MATRIX> > > flux_matrices_up;
00665 };
00666
00667
00668
00669
00670 namespace internal
00671 {
00672 template <class MATRIX>
00673 void
00674 reinit(MatrixBlock<MATRIX>& v, const BlockSparsityPattern& p)
00675 {
00676 v.row_indices = p.get_row_indices();
00677 v.column_indices = p.get_column_indices();
00678 }
00679
00680
00681 template <typename number>
00682 void
00683 reinit(MatrixBlock<::SparseMatrix<number> >& v, const BlockSparsityPattern& p)
00684 {
00685 v.row_indices = p.get_row_indices();
00686 v.column_indices = p.get_column_indices();
00687 v.matrix.reinit(p.block(v.row, v.column));
00688 }
00689 }
00690
00691
00692 template <class MATRIX>
00693 inline
00694 MatrixBlock<MATRIX>::MatrixBlock()
00695 :
00696 row(deal_II_numbers::invalid_unsigned_int),
00697 column(deal_II_numbers::invalid_unsigned_int)
00698 {}
00699
00700
00701 template <class MATRIX>
00702 inline
00703 MatrixBlock<MATRIX>::MatrixBlock(const MatrixBlock<MATRIX>& M)
00704 :
00705 Subscriptor(),
00706 row(M.row),
00707 column(M.column),
00708 matrix(M.matrix),
00709 row_indices(M.row_indices),
00710 column_indices(M.column_indices)
00711 {}
00712
00713
00714 template <class MATRIX>
00715 inline
00716 MatrixBlock<MATRIX>::MatrixBlock(unsigned int i, unsigned int j)
00717 :
00718 row(i), column(j)
00719 {}
00720
00721
00722 template <class MATRIX>
00723 inline
00724 void
00725 MatrixBlock<MATRIX>::reinit(const BlockSparsityPattern& sparsity)
00726 {
00727 internal::reinit(*this, sparsity);
00728 }
00729
00730
00731 template <class MATRIX>
00732 inline
00733 MatrixBlock<MATRIX>::operator MATRIX&()
00734 {
00735 return matrix;
00736 }
00737
00738
00739 template <class MATRIX>
00740 inline
00741 MatrixBlock<MATRIX>::operator const MATRIX&() const
00742 {
00743 return matrix;
00744 }
00745
00746
00747 template <class MATRIX>
00748 inline void
00749 MatrixBlock<MATRIX>::add (
00750 const unsigned int gi,
00751 const unsigned int gj,
00752 const typename MATRIX::value_type value)
00753 {
00754 Assert(row_indices.size() != 0, ExcNotInitialized());
00755 Assert(column_indices.size() != 0, ExcNotInitialized());
00756
00757 const std::pair<unsigned int, unsigned int> bi
00758 = row_indices.global_to_local(gi);
00759 const std::pair<unsigned int, unsigned int> bj
00760 = column_indices.global_to_local(gj);
00761
00762 Assert (bi.first == row, ExcBlockIndexMismatch(bi.first, row));
00763 Assert (bj.first == column, ExcBlockIndexMismatch(bj.first, column));
00764
00765 matrix.add(bi.second, bj.second, value);
00766 }
00767
00768
00769 template <class MATRIX>
00770 template <typename number>
00771 inline
00772 void
00773 MatrixBlock<MATRIX>::add (const std::vector<unsigned int>& r_indices,
00774 const std::vector<unsigned int>& c_indices,
00775 const FullMatrix<number>& values,
00776 const bool elide_zero_values)
00777 {
00778 Assert(row_indices.size() != 0, ExcNotInitialized());
00779 Assert(column_indices.size() != 0, ExcNotInitialized());
00780
00781 AssertDimension (r_indices.size(), values.m());
00782 AssertDimension (c_indices.size(), values.n());
00783
00784 for (unsigned int i=0; i<row_indices.size(); ++i)
00785 add (r_indices[i], c_indices.size(), &c_indices[0], &values(i,0),
00786 elide_zero_values);
00787 }
00788
00789
00790 template <class MATRIX>
00791 template <typename number>
00792 inline
00793 void
00794 MatrixBlock<MATRIX>::add (const unsigned int b_row,
00795 const unsigned int n_cols,
00796 const unsigned int *col_indices,
00797 const number *values,
00798 const bool,
00799 const bool)
00800 {
00801 Assert(row_indices.size() != 0, ExcNotInitialized());
00802 Assert(column_indices.size() != 0, ExcNotInitialized());
00803
00804 const std::pair<unsigned int, unsigned int> bi
00805 = row_indices.global_to_local(b_row);
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
00817
00818 for (unsigned int j=0;j<n_cols;++j)
00819 {
00820 const std::pair<unsigned int, unsigned int> bj
00821 = column_indices.global_to_local(col_indices[j]);
00822 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
00823
00824 matrix.add(bi.second, bj.second, values[j]);
00825 }
00826
00827 }
00828
00829
00830 template <class MATRIX>
00831 template <typename number>
00832 inline
00833 void
00834 MatrixBlock<MATRIX>::add (const std::vector<unsigned int> &indices,
00835 const FullMatrix<number> &values,
00836 const bool elide_zero_values)
00837 {
00838 Assert(row_indices.size() != 0, ExcNotInitialized());
00839 Assert(column_indices.size() != 0, ExcNotInitialized());
00840
00841 AssertDimension (indices.size(), values.m());
00842 Assert (values.n() == values.m(), ExcNotQuadratic());
00843
00844 for (unsigned int i=0; i<indices.size(); ++i)
00845 add (indices[i], indices.size(), &indices[0], &values(i,0),
00846 elide_zero_values);
00847 }
00848
00849
00850
00851 template <class MATRIX>
00852 template <typename number>
00853 inline
00854 void
00855 MatrixBlock<MATRIX>::add (const unsigned int row,
00856 const std::vector<unsigned int> &col_indices,
00857 const std::vector<number> &values,
00858 const bool elide_zero_values)
00859 {
00860 Assert(row_indices.size() != 0, ExcNotInitialized());
00861 Assert(column_indices.size() != 0, ExcNotInitialized());
00862
00863 AssertDimension (col_indices.size(), values.size());
00864 add (row, col_indices.size(), &col_indices[0], &values[0],
00865 elide_zero_values);
00866 }
00867
00868
00869 template <class MATRIX>
00870 template <class VECTOR>
00871 inline
00872 void
00873 MatrixBlock<MATRIX>::vmult (VECTOR& w, const VECTOR& v) const
00874 {
00875 matrix.vmult(w,v);
00876 }
00877
00878
00879 template <class MATRIX>
00880 template <class VECTOR>
00881 inline
00882 void
00883 MatrixBlock<MATRIX>::vmult_add (VECTOR& w, const VECTOR& v) const
00884 {
00885 matrix.vmult_add(w,v);
00886 }
00887
00888
00889 template <class MATRIX>
00890 template <class VECTOR>
00891 inline
00892 void
00893 MatrixBlock<MATRIX>::Tvmult (VECTOR& w, const VECTOR& v) const
00894 {
00895 matrix.Tvmult(w,v);
00896 }
00897
00898
00899 template <class MATRIX>
00900 template <class VECTOR>
00901 inline
00902 void
00903 MatrixBlock<MATRIX>::Tvmult_add (VECTOR& w, const VECTOR& v) const
00904 {
00905 matrix.Tvmult_add(w,v);
00906 }
00907
00908
00909 template <class MATRIX>
00910 inline
00911 std::size_t
00912 MatrixBlock<MATRIX>::memory_consumption () const
00913 {
00914 return (sizeof(*this)
00915 + MemoryConsumption::memory_consumption(matrix)
00916 - sizeof(matrix));
00917 }
00918
00919
00920
00921 template <class MATRIX>
00922 inline void
00923 MatrixBlockVector<MATRIX>::add(
00924 unsigned int row, unsigned int column,
00925 const std::string& name)
00926 {
00927 std_cxx1x::shared_ptr<value_type> p(new value_type(row, column));
00928 NamedData<std_cxx1x::shared_ptr<value_type> >::add(p, name);
00929 }
00930
00931
00932 template <class MATRIX>
00933 inline void
00934 MatrixBlockVector<MATRIX>::reinit(const BlockSparsityPattern& sparsity)
00935 {
00936 for (unsigned int i=0;i<this->size();++i)
00937 {
00938 block(i).reinit(sparsity);
00939 }
00940 }
00941
00942
00943 template <class MATRIX>
00944 inline void
00945 MatrixBlockVector<MATRIX>::clear(bool really_clean)
00946 {
00947 if (really_clean)
00948 {
00949 Assert(false, ExcNotImplemented());
00950 }
00951 else
00952 {
00953 for (unsigned int i=0;i<this->size();++i)
00954 matrix(i).clear();
00955 }
00956 }
00957
00958
00959
00960 template <class MATRIX>
00961 inline const MatrixBlock<MATRIX>&
00962 MatrixBlockVector<MATRIX>::block(unsigned int i) const
00963 {
00964 return *this->read(i);
00965 }
00966
00967
00968 template <class MATRIX>
00969 inline MatrixBlock<MATRIX>&
00970 MatrixBlockVector<MATRIX>::block(unsigned int i)
00971 {
00972 return *(*this)(i);
00973 }
00974
00975
00976 template <class MATRIX>
00977 inline MATRIX&
00978 MatrixBlockVector<MATRIX>::matrix(unsigned int i)
00979 {
00980 return (*this)(i)->matrix;
00981 }
00982
00983
00984
00985
00986
00987 template <class MATRIX>
00988 inline
00989 MGMatrixBlockVector<MATRIX>::MGMatrixBlockVector(
00990 const bool e, const bool f)
00991 :
00992 edge_matrices(e),
00993 edge_flux_matrices(f)
00994 {}
00995
00996
00997 template <class MATRIX>
00998 inline
00999 unsigned int
01000 MGMatrixBlockVector<MATRIX>::size () const
01001 {
01002 return matrices.size();
01003 }
01004
01005
01006 template <class MATRIX>
01007 inline void
01008 MGMatrixBlockVector<MATRIX>::add(
01009 unsigned int row, unsigned int column,
01010 const std::string& name)
01011 {
01012 MGLevelObject<MatrixBlock<MATRIX> > p(0, 1);
01013 p[0].row = row;
01014 p[0].column = column;
01015
01016 matrices.add(p, name);
01017 if (edge_matrices)
01018 {
01019 matrices_in.add(p, name);
01020 matrices_out.add(p, name);
01021 }
01022 if (edge_flux_matrices)
01023 {
01024 flux_matrices_up.add(p, name);
01025 flux_matrices_down.add(p, name);
01026 }
01027 }
01028
01029
01030 template <class MATRIX>
01031 inline const MGLevelObject<MatrixBlock<MATRIX> >&
01032 MGMatrixBlockVector<MATRIX>::block(unsigned int i) const
01033 {
01034 return matrices.read(i);
01035 }
01036
01037
01038 template <class MATRIX>
01039 inline MGLevelObject<MatrixBlock<MATRIX> >&
01040 MGMatrixBlockVector<MATRIX>::block(unsigned int i)
01041 {
01042 return matrices(i);
01043 }
01044
01045
01046 template <class MATRIX>
01047 inline const MGLevelObject<MatrixBlock<MATRIX> >&
01048 MGMatrixBlockVector<MATRIX>::block_in(unsigned int i) const
01049 {
01050 return matrices_in.read(i);
01051 }
01052
01053
01054 template <class MATRIX>
01055 inline MGLevelObject<MatrixBlock<MATRIX> >&
01056 MGMatrixBlockVector<MATRIX>::block_in(unsigned int i)
01057 {
01058 return matrices_in(i);
01059 }
01060
01061
01062 template <class MATRIX>
01063 inline const MGLevelObject<MatrixBlock<MATRIX> >&
01064 MGMatrixBlockVector<MATRIX>::block_out(unsigned int i) const
01065 {
01066 return matrices_out.read(i);
01067 }
01068
01069
01070 template <class MATRIX>
01071 inline MGLevelObject<MatrixBlock<MATRIX> >&
01072 MGMatrixBlockVector<MATRIX>::block_out(unsigned int i)
01073 {
01074 return matrices_out(i);
01075 }
01076
01077
01078 template <class MATRIX>
01079 inline const MGLevelObject<MatrixBlock<MATRIX> >&
01080 MGMatrixBlockVector<MATRIX>::block_up(unsigned int i) const
01081 {
01082 return flux_matrices_up.read(i);
01083 }
01084
01085
01086 template <class MATRIX>
01087 inline MGLevelObject<MatrixBlock<MATRIX> >&
01088 MGMatrixBlockVector<MATRIX>::block_up(unsigned int i)
01089 {
01090 return flux_matrices_up(i);
01091 }
01092
01093
01094 template <class MATRIX>
01095 inline const MGLevelObject<MatrixBlock<MATRIX> >&
01096 MGMatrixBlockVector<MATRIX>::block_down(unsigned int i) const
01097 {
01098 return flux_matrices_down.read(i);
01099 }
01100
01101
01102 template <class MATRIX>
01103 inline MGLevelObject<MatrixBlock<MATRIX> >&
01104 MGMatrixBlockVector<MATRIX>::block_down(unsigned int i)
01105 {
01106 return flux_matrices_down(i);
01107 }
01108
01109
01110 template <class MATRIX>
01111 inline void
01112 MGMatrixBlockVector<MATRIX>::reinit_matrix(const MGLevelObject<BlockSparsityPattern>& sparsity)
01113 {
01114 for (unsigned int i=0;i<this->size();++i)
01115 {
01116 MGLevelObject<MatrixBlock<MATRIX> >& o = block(i);
01117 const unsigned int row = o[o.min_level()].row;
01118 const unsigned int col = o[o.min_level()].column;
01119
01120 o.resize(sparsity.min_level(), sparsity.max_level());
01121 for (unsigned int level = o.min_level();level <= o.max_level();++level)
01122 {
01123 o[level].row = row;
01124 o[level].column = col;
01125 internal::reinit(o[level], sparsity[level]);
01126 }
01127 }
01128 }
01129
01130
01131 template <class MATRIX>
01132 inline void
01133 MGMatrixBlockVector<MATRIX>::reinit_edge(const MGLevelObject<BlockSparsityPattern>& sparsity)
01134 {
01135 for (unsigned int i=0;i<this->size();++i)
01136 {
01137 MGLevelObject<MatrixBlock<MATRIX> >& o = block(i);
01138 const unsigned int row = o[o.min_level()].row;
01139 const unsigned int col = o[o.min_level()].column;
01140
01141 block_in(i).resize(sparsity.min_level(), sparsity.max_level());
01142 block_out(i).resize(sparsity.min_level(), sparsity.max_level());
01143 for (unsigned int level = o.min_level();level <= o.max_level();++level)
01144 {
01145 block_in(i)[level].row = row;
01146 block_in(i)[level].column = col;
01147 internal::reinit(block_in(i)[level], sparsity[level]);
01148 block_out(i)[level].row = row;
01149 block_out(i)[level].column = col;
01150 internal::reinit(block_out(i)[level], sparsity[level]);
01151 }
01152 }
01153 }
01154
01155
01156 template <class MATRIX>
01157 inline void
01158 MGMatrixBlockVector<MATRIX>::reinit_edge_flux(const MGLevelObject<BlockSparsityPattern>& sparsity)
01159 {
01160 for (unsigned int i=0;i<this->size();++i)
01161 {
01162 MGLevelObject<MatrixBlock<MATRIX> >& o = block(i);
01163 const unsigned int row = o[o.min_level()].row;
01164 const unsigned int col = o[o.min_level()].column;
01165
01166 block_up(i).resize(sparsity.min_level(), sparsity.max_level());
01167 block_down(i).resize(sparsity.min_level(), sparsity.max_level());
01168 for (unsigned int level = o.min_level();level <= o.max_level();++level)
01169 {
01170 block_up(i)[level].row = row;
01171 block_up(i)[level].column = col;
01172 internal::reinit(block_up(i)[level], sparsity[level]);
01173 block_down(i)[level].row = row;
01174 block_down(i)[level].column = col;
01175 internal::reinit(block_down(i)[level], sparsity[level]);
01176 }
01177
01178 }
01179 }
01180
01181
01182 template <class MATRIX>
01183 inline void
01184 MGMatrixBlockVector<MATRIX>::clear_object(NamedData<MGLevelObject<MatrixBlock<MATRIX> > >& mo)
01185 {
01186 for (unsigned int i=0;i<mo.size();++i)
01187 {
01188 MGLevelObject<MatrixBlock<MATRIX> >& o = mo(i);
01189 for (unsigned int level = o.min_level();level <= o.max_level();++level)
01190 o[level].matrix.clear();
01191 }
01192 }
01193
01194
01195 template <class MATRIX>
01196 inline void
01197 MGMatrixBlockVector<MATRIX>::clear(bool really_clean)
01198 {
01199 if (really_clean)
01200 {
01201 Assert(false, ExcNotImplemented());
01202 }
01203 else
01204 {
01205 clear_object(matrices);
01206 clear_object(matrices_in);
01207 clear_object(matrices_out);
01208 clear_object(flux_matrices_up);
01209 clear_object(flux_matrices_down);
01210 }
01211 }
01212
01213
01214
01215
01216
01217 DEAL_II_NAMESPACE_CLOSE
01218
01219 #endif
01220