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