4#ifndef DUNE_ISTL_BVECTOR_HH
5#define DUNE_ISTL_BVECTOR_HH
10#include <initializer_list>
27#include "istlexception.hh"
46 template <
class B,
bool isNumber>
50 class BlockTraitsImp<B,true>
55 static constexpr unsigned int blockLevel()
62 class BlockTraitsImp<B,false>
65 using field_type =
typename B::field_type;
67 static constexpr unsigned int blockLevel()
76 using BlockTraits = BlockTraitsImp<B,IsNumber<B>::value>;
91 template<
class B,
class A=std::allocator<B> >
92 class block_vector_unmanaged :
public base_array_unmanaged<B,A>
97 using field_type =
typename Imp::BlockTraits<B>::field_type;
100 typedef B block_type;
103 typedef A allocator_type;
106 typedef typename A::size_type size_type;
109 typedef typename base_array_unmanaged<B,A>::iterator Iterator;
112 typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
115 typedef B value_type;
118 typedef B& reference;
121 typedef const B& const_reference;
126 block_vector_unmanaged& operator= (
const field_type& k)
128 for (size_type i=0; i<this->n; i++)
135 block_vector_unmanaged& operator+= (
const block_vector_unmanaged& y)
137#ifdef DUNE_ISTL_WITH_CHECKING
138 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
140 for (size_type i=0; i<this->n; ++i) (*
this)[i] += y[i];
145 block_vector_unmanaged& operator-= (
const block_vector_unmanaged& y)
147#ifdef DUNE_ISTL_WITH_CHECKING
148 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
150 for (size_type i=0; i<this->n; ++i) (*
this)[i] -= y[i];
155 block_vector_unmanaged& operator*= (
const field_type& k)
157 for (size_type i=0; i<this->n; ++i) (*
this)[i] *= k;
162 block_vector_unmanaged& operator/= (
const field_type& k)
164 for (size_type i=0; i<this->n; ++i) (*
this)[i] /= k;
169 block_vector_unmanaged& axpy (
const field_type& a,
const block_vector_unmanaged& y)
171#ifdef DUNE_ISTL_WITH_CHECKING
172 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
174 for (size_type i=0; i<this->n; ++i)
175 Impl::asVector((*
this)[i]).axpy(a,Impl::asVector(y[i]));
188 template<
class OtherB,
class OtherA>
189 auto operator* (
const block_vector_unmanaged<OtherB,OtherA>& y)
const
191 typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
193#ifdef DUNE_ISTL_WITH_CHECKING
194 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
196 for (size_type i=0; i<this->n; ++i) {
197 sum += PromotedType(((*
this)[i])*y[i]);
209 template<
class OtherB,
class OtherA>
210 auto dot(
const block_vector_unmanaged<OtherB,OtherA>& y)
const
212 typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
214#ifdef DUNE_ISTL_WITH_CHECKING
215 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
218 for (size_type i=0; i<this->n; ++i)
219 sum += Impl::asVector((*
this)[i]).dot(Impl::asVector(y[i]));
227 typename FieldTraits<field_type>::real_type one_norm ()
const
229 typename FieldTraits<field_type>::real_type sum=0;
230 for (size_type i=0; i<this->n; ++i)
231 sum += Impl::asVector((*
this)[i]).one_norm();
236 typename FieldTraits<field_type>::real_type one_norm_real ()
const
238 typename FieldTraits<field_type>::real_type sum=0;
239 for (size_type i=0; i<this->n; ++i)
240 sum += Impl::asVector((*
this)[i]).one_norm_real();
245 typename FieldTraits<field_type>::real_type two_norm ()
const
248 return sqrt(two_norm2());
252 typename FieldTraits<field_type>::real_type two_norm2 ()
const
254 typename FieldTraits<field_type>::real_type sum=0;
255 for (size_type i=0; i<this->n; ++i)
256 sum += Impl::asVector((*
this)[i]).two_norm2();
261 template <
typename ft = field_type,
262 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
263 typename FieldTraits<ft>::real_type infinity_norm()
const {
264 using real_type =
typename FieldTraits<ft>::real_type;
268 for (
auto const &xi : *
this) {
269 real_type
const a = Impl::asVector(xi).infinity_norm();
276 template <
typename ft = field_type,
277 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
278 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
279 using real_type =
typename FieldTraits<ft>::real_type;
283 for (
auto const &xi : *
this) {
284 real_type
const a = Impl::asVector(xi).infinity_norm_real();
291 template <
typename ft = field_type,
292 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
293 typename FieldTraits<ft>::real_type infinity_norm()
const {
294 using real_type =
typename FieldTraits<ft>::real_type;
301 for (
auto const &xi : *
this) {
302 real_type
const a = Impl::asVector(xi).infinity_norm();
306 return norm * (isNaN / isNaN);
310 template <
typename ft = field_type,
311 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
312 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
313 using real_type =
typename FieldTraits<ft>::real_type;
319 for (
auto const &xi : *
this) {
320 real_type
const a = Impl::asVector(xi).infinity_norm_real();
325 return norm * (isNaN / isNaN);
337 size_type dim ()
const
341 for (size_type i=0; i<this->n; i++)
342 d += Impl::asVector((*
this)[i]).dim();
349 block_vector_unmanaged () : base_array_unmanaged<B,A>()
363 ScopeGuard(F cleanupFunc) : cleanupFunc_(
std::move(cleanupFunc)) {}
364 ScopeGuard(
const ScopeGuard &) =
delete;
365 ScopeGuard(ScopeGuard &&) =
delete;
366 ScopeGuard &operator=(ScopeGuard) =
delete;
367 ~ScopeGuard() { cleanupFunc_(); }
381 ScopeGuard<F> makeScopeGuard(F cleanupFunc)
383 return { std::move(cleanupFunc) };
401 template<
class B,
class A=std::allocator<B> >
421 static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
424 typedef typename Imp::block_vector_unmanaged<B,A>::Iterator
Iterator;
427 typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator
ConstIterator;
464 static_assert(std::numeric_limits<S>::is_integer,
465 "capacity must be an unsigned integral type (be aware, that this constructor does not set the default value!)" );
467 storage_.reserve(_capacity);
485 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
509 "values are always copied")
523 return storage_.capacity();
539 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
540 storage_.resize(size);
561 "values are always copied")
572 noexcept(
noexcept(std::declval<BlockVector>().storage_ = a.storage_))
574 storage_ = a.storage_;
580 noexcept(
noexcept(std::declval<BlockVector>().swap(a)))
587 noexcept(
noexcept(std::declval<BlockVector>().storage_ = a.storage_))
590 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
591 storage_ = a.storage_;
597 noexcept(
noexcept(std::declval<BlockVector>().swap(a)))
606 std::declval<BlockVector&>().storage_.swap(other.storage_)))
608 DUNE_UNUSED const auto &guard = Imp::makeScopeGuard([&]{
610 other.syncBaseArray();
612 storage_.swap(other.storage_);
619 (
static_cast<Imp::block_vector_unmanaged<B,A>&
>(*this)) = k;
624 void syncBaseArray() noexcept
626 this->p = storage_.data();
627 this->n = storage_.size();
630 std::vector<B, A> storage_;
638 template<
class B,
class A>
639 struct FieldTraits< BlockVector<B, A> >
641 typedef typename FieldTraits<B>::field_type field_type;
642 typedef typename FieldTraits<B>::real_type real_type;
649 template<
class K,
class A>
654 for (size_type i=0; i<v.size(); i++)
655 s << v[i] << std::endl;
682 template<
class B,
class A>
684 template<
class B,
class A=std::allocator<B> >
686 class BlockVectorWindow :
public Imp::block_vector_unmanaged<B,A>
693 using field_type =
typename Imp::BlockTraits<B>::field_type;
696 typedef B block_type;
699 typedef A allocator_type;
702 typedef typename A::size_type size_type;
705 static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
708 typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
711 typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
716 BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
720 BlockVectorWindow (B* _p, size_type _n)
727 BlockVectorWindow (
const BlockVectorWindow& a)
734 BlockVectorWindow& operator= (
const BlockVectorWindow& a)
737#ifdef DUNE_ISTL_WITH_CHECKING
738 if (this->n!=a.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
744 for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
750 BlockVectorWindow& operator= (
const field_type& k)
752 (
static_cast<Imp::block_vector_unmanaged<B,A>&
>(*this)) = k;
757 operator BlockVector<B, A>()
const {
758 auto bv = BlockVector<B, A>(this->n);
760 std::copy(this->begin(), this->end(), bv.begin());
768 void set (size_type _n, B* _p)
775 void setsize (size_type _n)
793 size_type getsize ()
const
811 template<
class B,
class A=std::allocator<B> >
812 class compressed_block_vector_unmanaged :
public compressed_base_array_unmanaged<B,A>
819 using field_type =
typename Imp::BlockTraits<B>::field_type;
822 typedef B block_type;
825 typedef A allocator_type;
828 typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
831 typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
834 typedef typename A::size_type size_type;
838 compressed_block_vector_unmanaged& operator= (
const field_type& k)
840 for (size_type i=0; i<this->n; i++)
850 compressed_block_vector_unmanaged& operator+= (
const V& y)
852#ifdef DUNE_ISTL_WITH_CHECKING
853 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
855 for (size_type i=0; i<y.n; ++i) this->
operator[](y.j[i]) += y.p[i];
861 compressed_block_vector_unmanaged& operator-= (
const V& y)
863#ifdef DUNE_ISTL_WITH_CHECKING
864 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
866 for (size_type i=0; i<y.n; ++i) this->
operator[](y.j[i]) -= y.p[i];
872 compressed_block_vector_unmanaged& axpy (
const field_type& a,
const V& y)
874#ifdef DUNE_ISTL_WITH_CHECKING
875 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
877 for (size_type i=0; i<y.n; ++i)
878 Impl::asVector((*
this)[y.j[i]]).axpy(a,Impl::asVector(y.p[i]));
883 compressed_block_vector_unmanaged& operator*= (
const field_type& k)
885 for (size_type i=0; i<this->n; ++i) (this->p)[i] *= k;
890 compressed_block_vector_unmanaged& operator/= (
const field_type& k)
892 for (size_type i=0; i<this->n; ++i) (this->p)[i] /= k;
900 field_type operator* (
const compressed_block_vector_unmanaged& y)
const
902#ifdef DUNE_ISTL_WITH_CHECKING
903 if (!includesindexset(y) || !y.includesindexset(*
this) )
907 for (size_type i=0; i<this->n; ++i)
908 sum += (this->p)[i] * y[(this->j)[i]];
916 typename FieldTraits<field_type>::real_type one_norm ()
const
918 typename FieldTraits<field_type>::real_type sum=0;
919 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm();
924 typename FieldTraits<field_type>::real_type one_norm_real ()
const
926 typename FieldTraits<field_type>::real_type sum=0;
927 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm_real();
932 typename FieldTraits<field_type>::real_type two_norm ()
const
935 typename FieldTraits<field_type>::real_type sum=0;
936 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
941 typename FieldTraits<field_type>::real_type two_norm2 ()
const
943 typename FieldTraits<field_type>::real_type sum=0;
944 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
949 template <
typename ft = field_type,
950 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
951 typename FieldTraits<ft>::real_type infinity_norm()
const {
952 using real_type =
typename FieldTraits<ft>::real_type;
956 for (
auto const &x : *
this) {
957 real_type
const a = x.infinity_norm();
964 template <
typename ft = field_type,
965 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
966 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
967 using real_type =
typename FieldTraits<ft>::real_type;
971 for (
auto const &x : *
this) {
972 real_type
const a = x.infinity_norm_real();
979 template <
typename ft = field_type,
980 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
981 typename FieldTraits<ft>::real_type infinity_norm()
const {
982 using real_type =
typename FieldTraits<ft>::real_type;
987 for (
auto const &x : *
this) {
988 real_type
const a = x.infinity_norm();
992 return norm * (isNaN / isNaN);
996 template <
typename ft = field_type,
997 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
998 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
999 using real_type =
typename FieldTraits<ft>::real_type;
1003 real_type isNaN = 1;
1004 for (
auto const &x : *
this) {
1005 real_type
const a = x.infinity_norm_real();
1006 norm =
max(a, norm);
1009 return norm * (isNaN / isNaN);
1015 size_type N ()
const
1021 size_type dim ()
const
1024 for (size_type i=0; i<this->n; i++)
1025 d += (this->p)[i].dim();
1031 compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
1036 bool includesindexset (
const V& y)
1038 typename V::ConstIterator e=this->end();
1039 for (size_type i=0; i<y.n; i++)
1040 if (this->find(y.j[i])==e)
1065 template<
class B,
class A=std::allocator<B> >
1066 class CompressedBlockVectorWindow :
public compressed_block_vector_unmanaged<B,A>
1073 using field_type =
typename Imp::BlockTraits<B>::field_type;
1076 typedef B block_type;
1079 typedef A allocator_type;
1082 typedef typename A::size_type size_type;
1085 static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
1088 typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
1091 typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
1096 CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
1100 CompressedBlockVectorWindow (B* _p, size_type* _j, size_type _n)
1108 CompressedBlockVectorWindow (
const CompressedBlockVectorWindow& a)
1116 CompressedBlockVectorWindow& operator= (
const CompressedBlockVectorWindow& a)
1119#ifdef DUNE_ISTL_WITH_CHECKING
1120 if (this->n!=a.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
1126 for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
1127 for (size_type i=0; i<this->n; i++) this->j[i]=a.j[i];
1133 CompressedBlockVectorWindow& operator= (
const field_type& k)
1135 (
static_cast<compressed_block_vector_unmanaged<B,A>&
>(*this)) = k;
1143 void set (size_type _n, B* _p, size_type* _j)
1151 void setsize (size_type _n)
1163 void setindexptr (size_type* _j)
1175 size_type* getindexptr ()
1181 const B* getptr ()
const
1187 const size_type* getindexptr ()
const
1192 size_type getsize ()
const
1202 template<
typename B,
typename A>
1203 struct AutonomousValueType<Imp::BlockVectorWindow<B,A>>
1205 using type = BlockVector<B, A>;
Implements several basic array containers.
A vector of blocks with memory management.
Definition: bvector.hh:403
BlockVector()
makes empty vector
Definition: bvector.hh:432
void reserve(size_type capacity)
Reserve space.
Definition: bvector.hh:482
BlockVector(BlockVector &&a) noexcept(noexcept(std::declval< BlockVector >().swap(a)))
move constructor
Definition: bvector.hh:579
BlockVector(size_type _n)
make vector with _n components
Definition: bvector.hh:438
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:536
Imp::block_vector_unmanaged< B, A >::Iterator Iterator
make iterators available as types
Definition: bvector.hh:424
static constexpr unsigned int blocklevel
increment block level counter
Definition: bvector.hh:421
BlockVector(const BlockVector &a) noexcept(noexcept(std::declval< BlockVector >().storage_=a.storage_))
copy constructor
Definition: bvector.hh:571
A allocator_type
export the allocator type
Definition: bvector.hh:415
BlockVector & operator=(const BlockVector &a) noexcept(noexcept(std::declval< BlockVector >().storage_=a.storage_))
assignment
Definition: bvector.hh:586
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:409
A::size_type size_type
The type for the index access.
Definition: bvector.hh:418
BlockVector(std::initializer_list< B > const &l)
Construct from a std::initializer_list.
Definition: bvector.hh:444
size_type capacity() const
Get the capacity of the vector.
Definition: bvector.hh:521
BlockVector(size_type _n, S _capacity)
Make vector with _n components but preallocating capacity components.
Definition: bvector.hh:462
B block_type
export the type representing the components
Definition: bvector.hh:412
void swap(BlockVector &other) noexcept(noexcept(std::declval< BlockVector & >().storage_.swap(other.storage_)))
swap operation
Definition: bvector.hh:604
Imp::block_vector_unmanaged< B, A >::ConstIterator ConstIterator
make iterators available as types
Definition: bvector.hh:427
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
Provides the functions dot(a,b) := and dotT(a,b) := .
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Type traits to determine the type of reals (when working with complex numbers)
Implements a vector constructed from a given type representing a field and a compile-time given size.
auto dot(const A &a, const B &b) -> typename std::enable_if<!IsVector< A >::value &&!std::is_same< typename FieldTraits< A >::field_type, typename FieldTraits< A >::real_type > ::value, decltype(conj(a) *b)>::type
computes the dot product for fundamental data types according to Petsc's VectDot function: dot(a,...
Definition: dotproduct.hh:40
#define DUNE_UNUSED
A macro for marking variables that the compiler mistakenly flags as unused, which sometimes happens d...
Definition: unused.hh:16
#define DUNE_DEPRECATED_MSG(text)
Mark some entity as deprecated.
Definition: deprecated.hh:169
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:14
Implements a scalar vector view wrapper around an existing scalar.
Traits for type conversions and type information.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.