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>
57 class BlockTraitsImp<B,false>
60 using field_type =
typename B::field_type;
66 using BlockTraits = BlockTraitsImp<B,IsNumber<B>::value>;
81 template<
class B,
class A=std::allocator<B> >
82 class block_vector_unmanaged :
public base_array_unmanaged<B,A>
87 using field_type =
typename Imp::BlockTraits<B>::field_type;
93 typedef A allocator_type;
96 typedef typename A::size_type size_type;
99 typedef typename base_array_unmanaged<B,A>::iterator Iterator;
102 typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
105 typedef B value_type;
108 typedef B& reference;
111 typedef const B& const_reference;
116 block_vector_unmanaged& operator= (
const field_type& k)
118 for (size_type i=0; i<this->n; i++)
125 block_vector_unmanaged& operator+= (
const block_vector_unmanaged& y)
127#ifdef DUNE_ISTL_WITH_CHECKING
128 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
130 for (size_type i=0; i<this->n; ++i) (*
this)[i] += y[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 field_type& k)
147 for (size_type i=0; i<this->n; ++i) (*
this)[i] *= k;
152 block_vector_unmanaged& operator/= (
const field_type& k)
154 for (size_type i=0; i<this->n; ++i) (*
this)[i] /= k;
159 block_vector_unmanaged& axpy (
const field_type& a,
const block_vector_unmanaged& y)
161#ifdef DUNE_ISTL_WITH_CHECKING
162 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
164 for (size_type i=0; i<this->n; ++i)
165 Impl::asVector((*
this)[i]).axpy(a,Impl::asVector(y[i]));
178 template<
class OtherB,
class OtherA>
179 auto operator* (
const block_vector_unmanaged<OtherB,OtherA>& y)
const
181 typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
183#ifdef DUNE_ISTL_WITH_CHECKING
184 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
186 for (size_type i=0; i<this->n; ++i) {
187 sum += PromotedType(((*
this)[i])*y[i]);
199 template<
class OtherB,
class OtherA>
200 auto dot(
const block_vector_unmanaged<OtherB,OtherA>& y)
const
202 typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
204#ifdef DUNE_ISTL_WITH_CHECKING
205 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
208 for (size_type i=0; i<this->n; ++i)
209 sum += Impl::asVector((*
this)[i]).dot(Impl::asVector(y[i]));
217 typename FieldTraits<field_type>::real_type one_norm ()
const
219 typename FieldTraits<field_type>::real_type sum=0;
220 for (size_type i=0; i<this->n; ++i)
221 sum += Impl::asVector((*
this)[i]).one_norm();
226 typename FieldTraits<field_type>::real_type one_norm_real ()
const
228 typename FieldTraits<field_type>::real_type sum=0;
229 for (size_type i=0; i<this->n; ++i)
230 sum += Impl::asVector((*
this)[i]).one_norm_real();
235 typename FieldTraits<field_type>::real_type two_norm ()
const
238 return sqrt(two_norm2());
242 typename FieldTraits<field_type>::real_type two_norm2 ()
const
244 typename FieldTraits<field_type>::real_type sum=0;
245 for (size_type i=0; i<this->n; ++i)
246 sum += Impl::asVector((*
this)[i]).two_norm2();
251 template <
typename ft = field_type,
252 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
253 typename FieldTraits<ft>::real_type infinity_norm()
const {
254 using real_type =
typename FieldTraits<ft>::real_type;
258 for (
auto const &xi : *
this) {
259 real_type
const a = Impl::asVector(xi).infinity_norm();
266 template <
typename ft = field_type,
267 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
268 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
269 using real_type =
typename FieldTraits<ft>::real_type;
273 for (
auto const &xi : *
this) {
274 real_type
const a = Impl::asVector(xi).infinity_norm_real();
281 template <
typename ft = field_type,
282 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
283 typename FieldTraits<ft>::real_type infinity_norm()
const {
284 using real_type =
typename FieldTraits<ft>::real_type;
291 for (
auto const &xi : *
this) {
292 real_type
const a = Impl::asVector(xi).infinity_norm();
296 return norm * (isNaN / isNaN);
300 template <
typename ft = field_type,
301 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
302 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
303 using real_type =
typename FieldTraits<ft>::real_type;
309 for (
auto const &xi : *
this) {
310 real_type
const a = Impl::asVector(xi).infinity_norm_real();
315 return norm * (isNaN / isNaN);
327 size_type dim ()
const
331 for (size_type i=0; i<this->n; i++)
332 d += Impl::asVector((*
this)[i]).dim();
339 block_vector_unmanaged () : base_array_unmanaged<B,A>()
353 ScopeGuard(F cleanupFunc) : cleanupFunc_(
std::move(cleanupFunc)) {}
354 ScopeGuard(
const ScopeGuard &) =
delete;
355 ScopeGuard(ScopeGuard &&) =
delete;
356 ScopeGuard &operator=(ScopeGuard) =
delete;
357 ~ScopeGuard() { cleanupFunc_(); }
371 ScopeGuard<F> makeScopeGuard(F cleanupFunc)
373 return { std::move(cleanupFunc) };
391 template<
class B,
class A=std::allocator<B> >
411 [[deprecated(
"Use free function blockLevel(). Will be removed after 2.8.")]]
412 static constexpr unsigned int blocklevel = blockLevel<B>()+1;
415 typedef typename Imp::block_vector_unmanaged<B,A>::Iterator
Iterator;
418 typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator
ConstIterator;
455 static_assert(std::numeric_limits<S>::is_integer,
456 "capacity must be an unsigned integral type (be aware, that this constructor does not set the default value!)" );
458 storage_.reserve(_capacity);
475 [[maybe_unused]]
const auto &guard =
476 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
488 return storage_.capacity();
503 [[maybe_unused]]
const auto &guard =
504 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
505 storage_.resize(size);
510 noexcept(
noexcept(std::declval<BlockVector>().storage_ = a.storage_))
512 storage_ = a.storage_;
518 noexcept(
noexcept(std::declval<BlockVector>().swap(a)))
525 noexcept(
noexcept(std::declval<BlockVector>().storage_ = a.storage_))
527 [[maybe_unused]]
const auto &guard =
528 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
529 storage_ = a.storage_;
535 noexcept(
noexcept(std::declval<BlockVector>().swap(a)))
544 std::declval<BlockVector&>().storage_.swap(other.storage_)))
546 [[maybe_unused]]
const auto &guard = Imp::makeScopeGuard([&]{
548 other.syncBaseArray();
550 storage_.swap(other.storage_);
557 (
static_cast<Imp::block_vector_unmanaged<B,A>&
>(*this)) = k;
562 void syncBaseArray() noexcept
564 this->p = storage_.data();
565 this->n = storage_.size();
568 std::vector<B, A> storage_;
576 template<
class B,
class A>
577 struct FieldTraits< BlockVector<B, A> >
579 typedef typename FieldTraits<B>::field_type field_type;
580 typedef typename FieldTraits<B>::real_type real_type;
587 template<
class K,
class A>
592 for (size_type i=0; i<v.size(); i++)
593 s << v[i] << std::endl;
620 template<
class B,
class A>
622 template<
class B,
class A=std::allocator<B> >
624 class BlockVectorWindow :
public Imp::block_vector_unmanaged<B,A>
631 using field_type =
typename Imp::BlockTraits<B>::field_type;
634 typedef B block_type;
637 typedef A allocator_type;
640 typedef typename A::size_type size_type;
643 [[deprecated(
"Use free function blockLevel(). Will be removed after 2.8.")]]
644 static constexpr unsigned int blocklevel = blockLevel<B>()+1;
647 typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
650 typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
655 BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
659 BlockVectorWindow (B* _p, size_type _n)
666 BlockVectorWindow (
const BlockVectorWindow& a)
673 BlockVectorWindow& operator= (
const BlockVectorWindow& a)
676#ifdef DUNE_ISTL_WITH_CHECKING
677 if (this->n!=a.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
683 for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
689 BlockVectorWindow& operator= (
const field_type& k)
691 (
static_cast<Imp::block_vector_unmanaged<B,A>&
>(*this)) = k;
696 operator BlockVector<B, A>()
const {
697 auto bv = BlockVector<B, A>(this->n);
699 std::copy(this->begin(), this->end(), bv.begin());
707 void set (size_type _n, B* _p)
714 void setsize (size_type _n)
732 size_type getsize ()
const
750 template<
class B,
class A=std::allocator<B> >
751 class compressed_block_vector_unmanaged :
public compressed_base_array_unmanaged<B,A>
758 using field_type =
typename Imp::BlockTraits<B>::field_type;
761 typedef B block_type;
764 typedef A allocator_type;
767 typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
770 typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
773 typedef typename A::size_type size_type;
777 compressed_block_vector_unmanaged& operator= (
const field_type& k)
779 for (size_type i=0; i<this->n; i++)
789 compressed_block_vector_unmanaged& operator+= (
const V& y)
791#ifdef DUNE_ISTL_WITH_CHECKING
792 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
794 for (size_type i=0; i<y.n; ++i) this->
operator[](y.j[i]) += y.p[i];
800 compressed_block_vector_unmanaged& operator-= (
const V& y)
802#ifdef DUNE_ISTL_WITH_CHECKING
803 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
805 for (size_type i=0; i<y.n; ++i) this->
operator[](y.j[i]) -= y.p[i];
811 compressed_block_vector_unmanaged& axpy (
const field_type& a,
const V& y)
813#ifdef DUNE_ISTL_WITH_CHECKING
814 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
816 for (size_type i=0; i<y.n; ++i)
817 Impl::asVector((*
this)[y.j[i]]).axpy(a,Impl::asVector(y.p[i]));
822 compressed_block_vector_unmanaged& operator*= (
const field_type& k)
824 for (size_type i=0; i<this->n; ++i) (this->p)[i] *= k;
829 compressed_block_vector_unmanaged& operator/= (
const field_type& k)
831 for (size_type i=0; i<this->n; ++i) (this->p)[i] /= k;
839 field_type operator* (
const compressed_block_vector_unmanaged& y)
const
841#ifdef DUNE_ISTL_WITH_CHECKING
842 if (!includesindexset(y) || !y.includesindexset(*
this) )
846 for (size_type i=0; i<this->n; ++i)
847 sum += (this->p)[i] * y[(this->j)[i]];
855 typename FieldTraits<field_type>::real_type one_norm ()
const
857 typename FieldTraits<field_type>::real_type sum=0;
858 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm();
863 typename FieldTraits<field_type>::real_type one_norm_real ()
const
865 typename FieldTraits<field_type>::real_type sum=0;
866 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm_real();
871 typename FieldTraits<field_type>::real_type two_norm ()
const
874 typename FieldTraits<field_type>::real_type sum=0;
875 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
880 typename FieldTraits<field_type>::real_type two_norm2 ()
const
882 typename FieldTraits<field_type>::real_type sum=0;
883 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
888 template <
typename ft = field_type,
889 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
890 typename FieldTraits<ft>::real_type infinity_norm()
const {
891 using real_type =
typename FieldTraits<ft>::real_type;
895 for (
auto const &x : *
this) {
896 real_type
const a = x.infinity_norm();
903 template <
typename ft = field_type,
904 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
905 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
906 using real_type =
typename FieldTraits<ft>::real_type;
910 for (
auto const &x : *
this) {
911 real_type
const a = x.infinity_norm_real();
918 template <
typename ft = field_type,
919 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
920 typename FieldTraits<ft>::real_type infinity_norm()
const {
921 using real_type =
typename FieldTraits<ft>::real_type;
926 for (
auto const &x : *
this) {
927 real_type
const a = x.infinity_norm();
931 return norm * (isNaN / isNaN);
935 template <
typename ft = field_type,
936 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
937 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
938 using real_type =
typename FieldTraits<ft>::real_type;
943 for (
auto const &x : *
this) {
944 real_type
const a = x.infinity_norm_real();
948 return norm * (isNaN / isNaN);
960 size_type dim ()
const
963 for (size_type i=0; i<this->n; i++)
964 d += (this->p)[i].dim();
970 compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
975 bool includesindexset (
const V& y)
977 typename V::ConstIterator e=this->end();
978 for (size_type i=0; i<y.n; i++)
979 if (this->find(y.j[i])==e)
1004 template<
class B,
class A=std::allocator<B> >
1005 class CompressedBlockVectorWindow :
public compressed_block_vector_unmanaged<B,A>
1012 using field_type =
typename Imp::BlockTraits<B>::field_type;
1015 typedef B block_type;
1018 typedef A allocator_type;
1021 typedef typename A::size_type size_type;
1024 [[deprecated(
"Use free function blockLevel(). Will be removed after 2.8.")]]
1025 static constexpr unsigned int blocklevel = blockLevel<B>()+1;
1028 typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
1031 typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
1036 CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
1040 CompressedBlockVectorWindow (B* _p, size_type* _j, size_type _n)
1048 CompressedBlockVectorWindow (
const CompressedBlockVectorWindow& a)
1056 CompressedBlockVectorWindow& operator= (
const CompressedBlockVectorWindow& a)
1059#ifdef DUNE_ISTL_WITH_CHECKING
1060 if (this->n!=a.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
1066 for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
1067 for (size_type i=0; i<this->n; i++) this->j[i]=a.j[i];
1073 CompressedBlockVectorWindow& operator= (
const field_type& k)
1075 (
static_cast<compressed_block_vector_unmanaged<B,A>&
>(*this)) = k;
1083 void set (size_type _n, B* _p, size_type* _j)
1091 void setsize (size_type _n)
1103 void setindexptr (size_type* _j)
1115 size_type* getindexptr ()
1121 const B* getptr ()
const
1127 const size_type* getindexptr ()
const
1132 size_type getsize ()
const
1142 template<
typename B,
typename A>
1143 struct AutonomousValueType<Imp::BlockVectorWindow<B,A>>
1145 using type = BlockVector<B, A>;
Implements several basic array containers.
Helper functions for determining the vector/matrix block level.
A vector of blocks with memory management.
Definition: bvector.hh:393
BlockVector()
makes empty vector
Definition: bvector.hh:423
void reserve(size_type capacity)
Reserve space.
Definition: bvector.hh:473
BlockVector(BlockVector &&a) noexcept(noexcept(std::declval< BlockVector >().swap(a)))
move constructor
Definition: bvector.hh:517
BlockVector(size_type _n)
make vector with _n components
Definition: bvector.hh:429
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:501
Imp::block_vector_unmanaged< B, A >::Iterator Iterator
make iterators available as types
Definition: bvector.hh:415
static constexpr unsigned int blocklevel
increment block level counter
Definition: bvector.hh:412
BlockVector(const BlockVector &a) noexcept(noexcept(std::declval< BlockVector >().storage_=a.storage_))
copy constructor
Definition: bvector.hh:509
A allocator_type
export the allocator type
Definition: bvector.hh:405
BlockVector & operator=(const BlockVector &a) noexcept(noexcept(std::declval< BlockVector >().storage_=a.storage_))
assignment
Definition: bvector.hh:524
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:399
A::size_type size_type
The type for the index access.
Definition: bvector.hh:408
BlockVector(std::initializer_list< B > const &l)
Construct from a std::initializer_list.
Definition: bvector.hh:435
size_type capacity() const
Get the capacity of the vector.
Definition: bvector.hh:486
BlockVector(size_type _n, S _capacity)
Make vector with _n components but preallocating capacity components.
Definition: bvector.hh:453
B block_type
export the type representing the components
Definition: bvector.hh:402
void swap(BlockVector &other) noexcept(noexcept(std::declval< BlockVector & >().storage_.swap(other.storage_)))
swap operation
Definition: bvector.hh:542
Imp::block_vector_unmanaged< B, A >::ConstIterator ConstIterator
make iterators available as types
Definition: bvector.hh:418
Provides the functions dot(a,b) := and dotT(a,b) := .
Traits for type conversions and type information.
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_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:11
Implements a scalar vector view wrapper around an existing scalar.