1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_BACKEND_ISTL_BCRSMATRIX_HH
3 #define DUNE_PDELAB_BACKEND_ISTL_BCRSMATRIX_HH
4 
5 #include <dune/common/typetraits.hh>
6 #include <dune/common/shared_ptr.hh>
7 #if DUNE_VERSION_GTE(ISTL,2,8)
8 #include <dune/istl/blocklevel.hh>
9 #endif
10 #include <dune/pdelab/backend/common/tags.hh>
11 #include <dune/pdelab/backend/common/uncachedmatrixview.hh>
12 #include <dune/pdelab/backend/common/aliasedmatrixview.hh>
13 #include <dune/pdelab/backend/istl/matrixhelpers.hh>
14 #include <dune/pdelab/backend/istl/descriptors.hh>
15 
16 namespace Dune {
17   namespace PDELab {
18 
19     namespace ISTL {
20 
21       template<typename GFSV, typename GFSU, typename C, typename Stats>
22       class BCRSMatrix
23         : public Backend::impl::Wrapper<C>
24       {
25 
26         friend Backend::impl::Wrapper<C>;
27 
28       public:
29 
30         typedef typename C::field_type ElementType;
31         typedef ElementType E;
32         typedef C Container;
33         typedef typename C::field_type field_type;
34         typedef typename C::block_type block_type;
35         typedef typename C::size_type size_type;
36 
37         typedef GFSU TrialGridFunctionSpace;
38         typedef GFSV TestGridFunctionSpace;
39 
40         typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
41         typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
42 
43         typedef typename ISTL::build_pattern_type<C,GFSV,GFSU,typename GFSV::Ordering::ContainerAllocationTag>::type Pattern;
44 
45         typedef Stats PatternStatistics;
46 
47         using value_type = E;
48 
49 #ifndef DOXYGEN
50 
51         // some trickery to avoid exposing average users to the fact that there might
52         // be multiple statistics objects
53         typedef typename std::conditional<
54 #if DUNE_VERSION_LT(DUNE_ISTL,2,8)
55           (C::blocklevel > 2),
56 #else
57           (blockLevel<C>() > 2),
58 #endif
59           std::vector<PatternStatistics>,
60           PatternStatistics
61           >::type StatisticsReturnType;
62 
63 #endif // DOXYGEN
64 
65         template<typename RowCache, typename ColCache>
66         using LocalView = UncachedMatrixView<BCRSMatrix,RowCache,ColCache>;
67 
68         template<typename RowCache, typename ColCache>
69         using ConstLocalView = ConstUncachedMatrixView<const BCRSMatrix,RowCache,ColCache>;
70 
71         template<typename RowCache, typename ColCache>
72         using AliasedLocalView = AliasedMatrixView<BCRSMatrix,RowCache,ColCache>;
73 
74         template<typename RowCache, typename ColCache>
75         using ConstAliasedLocalView = ConstAliasedMatrixView<const BCRSMatrix,RowCache,ColCache>;
76 
77         template<typename GO>
BCRSMatrix(const GO & go)78         explicit BCRSMatrix (const GO& go)
79           : _container(std::make_shared<Container>())
80         {
81           _stats = go.matrixBackend().buildPattern(go,*this);
82         }
83 
84         /** \brief Construct matrix container using an externally given matrix as storage
85          *
86          * \tparam GO GridOperator type used to assemble into the matrix
87          *
88          * \param go GridOperator object used to assemble into the matrix
89          * \param container ISTL matrix type that stores the actual data
90          *
91          * This BCRSMatrix constructor will reassemble the matrix occupation pattern.
92          */
93         template<typename GO>
BCRSMatrix(const GO & go,Container & container)94         BCRSMatrix (const GO& go, Container& container)
95           : _container(Dune::stackobject_to_shared_ptr(container))
96         {
97           _stats = go.matrixBackend().buildPattern(go,*this);
98         }
99 
100         template<typename GO>
BCRSMatrix(const GO & go,const E & e)101         BCRSMatrix (const GO& go, const E& e)
102           : _container(std::make_shared<Container>())
103         {
104           _stats = go.matrixBackend().buildPattern(go,*this);
105           (*_container) = e;
106         }
107 
108         //! Creates an BCRSMatrix without allocating an underlying ISTL matrix.
BCRSMatrix(Backend::unattached_container=Backend::unattached_container ())109         explicit BCRSMatrix (Backend::unattached_container = Backend::unattached_container())
110         {}
111 
112         //! Creates an BCRSMatrix with an empty underlying ISTL matrix.
BCRSMatrix(Backend::attached_container)113         explicit BCRSMatrix (Backend::attached_container)
114           : _container(std::make_shared<Container>())
115         {}
116 
BCRSMatrix(const BCRSMatrix & rhs)117         BCRSMatrix(const BCRSMatrix& rhs)
118           : _container(std::make_shared<Container>(*(rhs._container)))
119         {}
120 
operator =(const BCRSMatrix & rhs)121         BCRSMatrix& operator=(const BCRSMatrix& rhs)
122         {
123           if (this == &rhs)
124             return *this;
125           _stats.clear();
126           if (attached())
127             {
128               (*_container) = (*(rhs._container));
129             }
130           else
131             {
132               _container = std::make_shared<Container>(*(rhs._container));
133             }
134           return *this;
135         }
136 
137         //! Returns pattern statistics for all contained BCRSMatrix objects.
patternStatistics() const138         const StatisticsReturnType& patternStatistics() const
139         {
140 #if DUNE_VERSION_LT(DUNE_ISTL,2,8)
141           return patternStatistics(std::integral_constant<bool,(C::blocklevel > 2)>());
142 #else
143           return patternStatistics(std::integral_constant<bool,(blockLevel<C>() > 2)>());
144 #endif
145         }
146 
147 #ifndef DOXYGEN
148 
149       private:
150 
patternStatistics(std::false_type multiple) const151         const PatternStatistics& patternStatistics(std::false_type multiple) const
152         {
153           if (_stats.empty())
154             DUNE_THROW(InvalidStateException,"no pattern statistics available");
155           return _stats[0];
156         }
157 
patternStatistics(std::true_type multiple) const158         const std::vector<PatternStatistics>& patternStatistics(std::true_type multiple) const
159         {
160           if (_stats.empty())
161             DUNE_THROW(InvalidStateException,"no pattern statistics available");
162           return _stats;
163         }
164 
165       public:
166 
167 #endif
168 
detach()169         void detach()
170         {
171           _container.reset();
172           _stats.clear();
173         }
174 
attach(std::shared_ptr<Container> container)175         void attach(std::shared_ptr<Container> container)
176         {
177           _container = container;
178         }
179 
attached() const180         bool attached() const
181         {
182           return bool(_container);
183         }
184 
storage() const185         const std::shared_ptr<Container>& storage() const
186         {
187           return _container;
188         }
189 
N() const190         size_type N() const
191         {
192           return _container->N();
193         }
194 
M() const195         size_type M() const
196         {
197           return _container->M();
198         }
199 
operator =(const E & e)200         BCRSMatrix& operator= (const E& e)
201         {
202           (*_container) = e;
203           return *this;
204         }
205 
operator *=(const E & e)206         BCRSMatrix& operator*= (const E& e)
207         {
208           (*_container) *= e;
209           return *this;
210         }
211 
operator ()(const RowIndex & ri,const ColIndex & ci)212         E& operator()(const RowIndex& ri, const ColIndex& ci)
213         {
214           return ISTL::access_matrix_element(ISTL::container_tag(*_container),*_container,ri,ci,ri.size()-1,ci.size()-1);
215         }
216 
operator ()(const RowIndex & ri,const ColIndex & ci) const217         const E& operator()(const RowIndex& ri, const ColIndex& ci) const
218         {
219           return ISTL::access_matrix_element(ISTL::container_tag(*_container),*_container,ri,ci,ri.size()-1,ci.size()-1);
220         }
221 
222       private:
223 
native() const224         const Container& native() const
225         {
226           return *_container;
227         }
228 
native()229         Container& native()
230         {
231           return *_container;
232         }
233 
234       public:
235 
236         template<typename RowCache, typename ColCache>
data(const RowCache & row_cache,const ColCache & col_cache)237         value_type* data(const RowCache& row_cache, const ColCache& col_cache)
238         {
239           return &((*this)(row_cache.containerIndex(0),col_cache.containerIndex(0)));
240         }
241 
242         template<typename RowCache, typename ColCache>
data(const RowCache & row_cache,const ColCache & col_cache) const243         const value_type* data(const RowCache& row_cache, const ColCache& col_cache) const
244         {
245           return &((*this)(row_cache.containerIndex(0),col_cache.containerIndex(0)));
246         }
247 
flush()248         void flush()
249         {}
250 
finalize()251         void finalize()
252         {}
253 
clear_row(const RowIndex & ri,const E & diagonal_entry)254         void clear_row(const RowIndex& ri, const E& diagonal_entry)
255         {
256           ISTL::clear_matrix_row(ISTL::container_tag(*_container),*_container,ri,ri.size()-1);
257           ISTL::write_matrix_element_if_exists(diagonal_entry,ISTL::container_tag(*_container),*_container,ri,ri,ri.size()-1,ri.size()-1);
258         }
259 
clear_row_block(const RowIndex & ri,const E & diagonal_entry)260         void clear_row_block(const RowIndex& ri, const E& diagonal_entry)
261         {
262           ISTL::clear_matrix_row_block(ISTL::container_tag(*_container),*_container,ri,ri.size()-1);
263           ISTL::write_matrix_element_if_exists_to_block(diagonal_entry,ISTL::container_tag(*_container),*_container,ri,ri,ri.size()-1,ri.size()-1);
264         }
265 
266       private:
267 
268         std::shared_ptr<Container> _container;
269         std::vector<PatternStatistics> _stats;
270 
271       };
272 
273     } // namespace ISTL
274 
275   } // namespace PDELab
276 } // namespace Dune
277 
278 #endif // DUNE_PDELAB_BACKEND_ISTL_BCRSMATRIX_HH
279