6#ifndef DUNE_ISTL_BVECTOR_HH
7#define DUNE_ISTL_BVECTOR_HH
12#include <initializer_list>
29#include "istlexception.hh"
48 template <
class B,
bool isNumber>
52 class BlockTraitsImp<B,true>
59 class BlockTraitsImp<B,false>
62 using field_type =
typename B::field_type;
68 using BlockTraits = BlockTraitsImp<B,IsNumber<B>::value>;
83 template<
class B,
class ST=std::
size_t >
84 class block_vector_unmanaged :
public base_array_unmanaged<B,ST>
89 using field_type =
typename Imp::BlockTraits<B>::field_type;
98 typedef typename base_array_unmanaged<B,ST>::iterator Iterator;
101 typedef typename base_array_unmanaged<B,ST>::const_iterator ConstIterator;
104 typedef B value_type;
107 typedef B& reference;
110 typedef const B& const_reference;
115 block_vector_unmanaged& operator= (
const field_type& k)
117 for (size_type i=0; i<this->n; i++)
124 block_vector_unmanaged& operator+= (
const block_vector_unmanaged& y)
126#ifdef DUNE_ISTL_WITH_CHECKING
127 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
129 for (size_type i=0; i<this->n; ++i) (*
this)[i] += y[i];
134 block_vector_unmanaged& operator-= (
const block_vector_unmanaged& y)
136#ifdef DUNE_ISTL_WITH_CHECKING
137 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
139 for (size_type i=0; i<this->n; ++i) (*
this)[i] -= y[i];
144 block_vector_unmanaged& operator*= (
const field_type& k)
146 for (size_type i=0; i<this->n; ++i) (*
this)[i] *= k;
151 block_vector_unmanaged& operator/= (
const field_type& k)
153 for (size_type i=0; i<this->n; ++i) (*
this)[i] /= k;
158 block_vector_unmanaged& axpy (
const field_type& a,
const block_vector_unmanaged& y)
160#ifdef DUNE_ISTL_WITH_CHECKING
161 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
163 for (size_type i=0; i<this->n; ++i)
164 Impl::asVector((*
this)[i]).axpy(a,Impl::asVector(y[i]));
177 template<
class OtherB,
class OtherST>
178 auto operator* (
const block_vector_unmanaged<OtherB,OtherST>& y)
const
180 typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
182#ifdef DUNE_ISTL_WITH_CHECKING
183 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
185 for (size_type i=0; i<this->n; ++i) {
186 sum += PromotedType(((*
this)[i])*y[i]);
198 template<
class OtherB,
class OtherST>
199 auto dot(
const block_vector_unmanaged<OtherB,OtherST>& y)
const
201 typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
203#ifdef DUNE_ISTL_WITH_CHECKING
204 if (this->n!=y.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
207 for (size_type i=0; i<this->n; ++i)
208 sum += Impl::asVector((*
this)[i]).dot(Impl::asVector(y[i]));
216 typename FieldTraits<field_type>::real_type one_norm ()
const
218 typename FieldTraits<field_type>::real_type sum=0;
219 for (size_type i=0; i<this->n; ++i)
220 sum += Impl::asVector((*
this)[i]).one_norm();
225 typename FieldTraits<field_type>::real_type one_norm_real ()
const
227 typename FieldTraits<field_type>::real_type sum=0;
228 for (size_type i=0; i<this->n; ++i)
229 sum += Impl::asVector((*
this)[i]).one_norm_real();
234 typename FieldTraits<field_type>::real_type two_norm ()
const
237 return sqrt(two_norm2());
241 typename FieldTraits<field_type>::real_type two_norm2 ()
const
243 typename FieldTraits<field_type>::real_type sum=0;
244 for (size_type i=0; i<this->n; ++i)
245 sum += Impl::asVector((*
this)[i]).two_norm2();
250 template <
typename ft = field_type,
251 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
252 typename FieldTraits<ft>::real_type infinity_norm()
const {
253 using real_type =
typename FieldTraits<ft>::real_type;
257 for (
auto const &xi : *
this) {
258 real_type
const a = Impl::asVector(xi).infinity_norm();
265 template <
typename ft = field_type,
266 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
267 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
268 using real_type =
typename FieldTraits<ft>::real_type;
272 for (
auto const &xi : *
this) {
273 real_type
const a = Impl::asVector(xi).infinity_norm_real();
280 template <
typename ft = field_type,
281 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
282 typename FieldTraits<ft>::real_type infinity_norm()
const {
283 using real_type =
typename FieldTraits<ft>::real_type;
290 for (
auto const &xi : *
this) {
291 real_type
const a = Impl::asVector(xi).infinity_norm();
299 template <
typename ft = field_type,
300 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
301 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
302 using real_type =
typename FieldTraits<ft>::real_type;
308 for (
auto const &xi : *
this) {
309 real_type
const a = Impl::asVector(xi).infinity_norm_real();
326 size_type dim ()
const
330 for (size_type i=0; i<this->n; i++)
331 d += Impl::asVector((*
this)[i]).dim();
338 block_vector_unmanaged () : base_array_unmanaged<B,ST>()
352 ScopeGuard(F cleanupFunc) : cleanupFunc_(
std::move(cleanupFunc)) {}
353 ScopeGuard(
const ScopeGuard &) =
delete;
354 ScopeGuard(ScopeGuard &&) =
delete;
355 ScopeGuard &operator=(ScopeGuard) =
delete;
356 ~ScopeGuard() { cleanupFunc_(); }
370 ScopeGuard<F> makeScopeGuard(F cleanupFunc)
372 return { std::move(cleanupFunc) };
390 template<
class B,
class A=std::allocator<B> >
391 class BlockVector :
public Imp::block_vector_unmanaged<B,typename std::allocator_traits<A>::size_type>
407 typedef typename std::allocator_traits<A>::size_type
size_type;
410 typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator
Iterator;
413 typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator
ConstIterator;
450 static_assert(std::numeric_limits<S>::is_integer,
451 "capacity must be an unsigned integral type (be aware, that this constructor does not set the default value!)" );
453 storage_.reserve(_capacity);
470 [[maybe_unused]]
const auto &guard =
471 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
483 return storage_.capacity();
498 [[maybe_unused]]
const auto &guard =
499 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
500 storage_.resize(
size);
505 noexcept(
noexcept(std::declval<BlockVector>().storage_ = a.storage_))
507 storage_ = a.storage_;
513 noexcept(
noexcept(std::declval<BlockVector>().swap(a)))
520 noexcept(
noexcept(std::declval<BlockVector>().storage_ = a.storage_))
522 [[maybe_unused]]
const auto &guard =
523 Imp::makeScopeGuard([
this]{ syncBaseArray(); });
524 storage_ = a.storage_;
530 noexcept(
noexcept(std::declval<BlockVector>().swap(a)))
539 std::declval<BlockVector&>().storage_.swap(other.storage_)))
541 [[maybe_unused]]
const auto &guard = Imp::makeScopeGuard([&]{
543 other.syncBaseArray();
545 storage_.swap(other.storage_);
552 (
static_cast<Imp::block_vector_unmanaged<B,size_type>&
>(*this)) = k;
557 void syncBaseArray() noexcept
559 this->p = storage_.data();
560 this->n = storage_.size();
563 using block_allocator_t =
typename std::allocator_traits<A>::template rebind_alloc<B>;
565 std::vector<B, block_allocator_t> storage_;
573 template<
class B,
class A>
574 struct FieldTraits< BlockVector<B, A> >
576 typedef typename FieldTraits<B>::field_type field_type;
577 typedef typename FieldTraits<B>::real_type real_type;
584 template<
class K,
class A>
589 for (size_type i=0; i<v.size(); i++)
590 s << v[i] << std::endl;
617 template<
class B,
class A>
619 template<
class B,
class A=std::allocator<B> >
621 class BlockVectorWindow :
public Imp::block_vector_unmanaged<B,typename A::size_type>
628 using field_type =
typename Imp::BlockTraits<B>::field_type;
631 typedef B block_type;
634 typedef A allocator_type;
637 typedef typename A::size_type size_type;
640 typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
643 typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
648 BlockVectorWindow () : Imp::block_vector_unmanaged<B,size_type>()
652 BlockVectorWindow (B* _p, size_type _n)
659 BlockVectorWindow (
const BlockVectorWindow& a)
666 BlockVectorWindow& operator= (
const BlockVectorWindow& a)
669#ifdef DUNE_ISTL_WITH_CHECKING
670 if (this->n!=a.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
676 for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
682 BlockVectorWindow& operator= (
const field_type& k)
684 (
static_cast<Imp::block_vector_unmanaged<B,size_type>&
>(*this)) = k;
689 operator BlockVector<B, A>()
const {
690 auto bv = BlockVector<B, A>(this->n);
692 std::copy(this->begin(), this->end(), bv.begin());
700 void set (size_type _n, B* _p)
707 void setsize (size_type _n)
725 size_type getsize ()
const
743 template<
class B,
class ST=std::
size_t >
744 class compressed_block_vector_unmanaged :
public compressed_base_array_unmanaged<B,ST>
751 using field_type =
typename Imp::BlockTraits<B>::field_type;
754 typedef B block_type;
757 typedef typename compressed_base_array_unmanaged<B,ST>::iterator Iterator;
760 typedef typename compressed_base_array_unmanaged<B,ST>::const_iterator ConstIterator;
763 typedef ST size_type;
767 compressed_block_vector_unmanaged& operator= (
const field_type& k)
769 for (size_type i=0; i<this->n; i++)
779 compressed_block_vector_unmanaged& operator+= (
const V& y)
781#ifdef DUNE_ISTL_WITH_CHECKING
782 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
784 for (size_type i=0; i<y.n; ++i) this->
operator[](y.j[i]) += y.p[i];
790 compressed_block_vector_unmanaged& operator-= (
const V& y)
792#ifdef DUNE_ISTL_WITH_CHECKING
793 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
795 for (size_type i=0; i<y.n; ++i) this->
operator[](y.j[i]) -= y.p[i];
801 compressed_block_vector_unmanaged& axpy (
const field_type& a,
const V& y)
803#ifdef DUNE_ISTL_WITH_CHECKING
804 if (!includesindexset(y))
DUNE_THROW(ISTLError,
"index set mismatch");
806 for (size_type i=0; i<y.n; ++i)
807 Impl::asVector((*
this)[y.j[i]]).axpy(a,Impl::asVector(y.p[i]));
812 compressed_block_vector_unmanaged& operator*= (
const field_type& k)
814 for (size_type i=0; i<this->n; ++i) (this->p)[i] *= k;
819 compressed_block_vector_unmanaged& operator/= (
const field_type& k)
821 for (size_type i=0; i<this->n; ++i) (this->p)[i] /= k;
829 field_type operator* (
const compressed_block_vector_unmanaged& y)
const
831#ifdef DUNE_ISTL_WITH_CHECKING
832 if (!includesindexset(y) || !y.includesindexset(*
this) )
836 for (size_type i=0; i<this->n; ++i)
837 sum += (this->p)[i] * y[(this->j)[i]];
845 typename FieldTraits<field_type>::real_type one_norm ()
const
847 typename FieldTraits<field_type>::real_type sum=0;
848 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm();
853 typename FieldTraits<field_type>::real_type one_norm_real ()
const
855 typename FieldTraits<field_type>::real_type sum=0;
856 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm_real();
861 typename FieldTraits<field_type>::real_type two_norm ()
const
864 typename FieldTraits<field_type>::real_type sum=0;
865 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
870 typename FieldTraits<field_type>::real_type two_norm2 ()
const
872 typename FieldTraits<field_type>::real_type sum=0;
873 for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
878 template <
typename ft = field_type,
879 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
880 typename FieldTraits<ft>::real_type infinity_norm()
const {
881 using real_type =
typename FieldTraits<ft>::real_type;
885 for (
auto const &x : *
this) {
886 real_type
const a = x.infinity_norm();
893 template <
typename ft = field_type,
894 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
895 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
896 using real_type =
typename FieldTraits<ft>::real_type;
900 for (
auto const &x : *
this) {
901 real_type
const a = x.infinity_norm_real();
908 template <
typename ft = field_type,
909 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
910 typename FieldTraits<ft>::real_type infinity_norm()
const {
911 using real_type =
typename FieldTraits<ft>::real_type;
916 for (
auto const &x : *
this) {
917 real_type
const a = x.infinity_norm();
925 template <
typename ft = field_type,
926 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
927 typename FieldTraits<ft>::real_type infinity_norm_real()
const {
928 using real_type =
typename FieldTraits<ft>::real_type;
933 for (
auto const &x : *
this) {
934 real_type
const a = x.infinity_norm_real();
950 size_type dim ()
const
953 for (size_type i=0; i<this->n; i++)
954 d += (this->p)[i].dim();
960 compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,ST>()
965 bool includesindexset (
const V& y)
967 typename V::ConstIterator e=this->end();
968 for (size_type i=0; i<y.n; i++)
969 if (this->find(y.j[i])==e)
997 template<
class B,
class ST=std::
size_t >
998 class CompressedBlockVectorWindow :
public compressed_block_vector_unmanaged<B,ST>
1005 using field_type =
typename Imp::BlockTraits<B>::field_type;
1008 typedef B block_type;
1011 typedef ST size_type;
1014 typedef typename compressed_block_vector_unmanaged<B,ST>::Iterator Iterator;
1017 typedef typename compressed_block_vector_unmanaged<B,ST>::ConstIterator ConstIterator;
1022 CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,ST>()
1026 CompressedBlockVectorWindow (B* _p, size_type* _j, size_type _n)
1033 [[deprecated(
"Reference semantics (i.e., assignment does shallow copy) will be removed due to ambiguity. "
1034 "Instead, use other explicit constructors or populate objects with the setter methods. "
1035 "This will be removed after Dune 2.11.")]]
1036 CompressedBlockVectorWindow(
const CompressedBlockVectorWindow &a) =
default;
1038 [[deprecated(
"Reference semantics (i.e., assignment does shallow copy) will be removed due to ambiguity. "
1039 "Instead, use other explicit constructors or populate objects with the setter methods. "
1040 "This will be removed after Dune 2.11.")]]
1041 CompressedBlockVectorWindow (CompressedBlockVectorWindow&& a) =
default;
1044 CompressedBlockVectorWindow& operator= (
const CompressedBlockVectorWindow& a)
1047#ifdef DUNE_ISTL_WITH_CHECKING
1048 if (this->n!=a.N())
DUNE_THROW(ISTLError,
"vector size mismatch");
1054 for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
1055 for (size_type i=0; i<this->n; i++) this->j[i]=a.j[i];
1061 CompressedBlockVectorWindow& operator= (
const field_type& k)
1063 (
static_cast<compressed_block_vector_unmanaged<B,ST>&
>(*this)) = k;
1067 [[deprecated(
"Reference semantics (i.e., assignment does shallow copy) will be removed due to ambiguity. "
1068 "Instead, use other explicit constructors or populate objects with the setter methods. "
1069 "This will be removed after Dune 2.11.")]]
1070 CompressedBlockVectorWindow& operator= (CompressedBlockVectorWindow&& a)
1072 return (*
this = std::as_const(a));
1078 void set (size_type _n, B* _p, size_type* _j)
1086 void setsize (size_type _n)
1098 void setindexptr (size_type* _j)
1110 size_type* getindexptr ()
1116 const B* getptr ()
const
1122 const size_type* getindexptr ()
const
1127 size_type getsize ()
const
1137 template<
typename B,
typename A>
1138 struct AutonomousValueType<Imp::BlockVectorWindow<B,A>>
1140 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:392
BlockVector()
makes empty vector
Definition: bvector.hh:418
Imp::block_vector_unmanaged< B, size_type >::ConstIterator ConstIterator
make iterators available as types
Definition: bvector.hh:413
void reserve(size_type capacity)
Reserve space.
Definition: bvector.hh:468
BlockVector(BlockVector &&a) noexcept(noexcept(std::declval< BlockVector >().swap(a)))
move constructor
Definition: bvector.hh:512
BlockVector(size_type _n)
make vector with _n components
Definition: bvector.hh:424
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:496
Imp::block_vector_unmanaged< B, size_type >::Iterator Iterator
make iterators available as types
Definition: bvector.hh:410
BlockVector(const BlockVector &a) noexcept(noexcept(std::declval< BlockVector >().storage_=a.storage_))
copy constructor
Definition: bvector.hh:504
A allocator_type
export the allocator type
Definition: bvector.hh:404
BlockVector & operator=(const BlockVector &a) noexcept(noexcept(std::declval< BlockVector >().storage_=a.storage_))
assignment
Definition: bvector.hh:519
std::allocator_traits< A >::size_type size_type
The type for the index access.
Definition: bvector.hh:407
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:398
BlockVector(std::initializer_list< B > const &l)
Construct from a std::initializer_list.
Definition: bvector.hh:430
size_type capacity() const
Get the capacity of the vector.
Definition: bvector.hh:481
BlockVector(size_type _n, S _capacity)
Make vector with _n components but preallocating capacity components.
Definition: bvector.hh:448
B block_type
export the type representing the components
Definition: bvector.hh:401
void swap(BlockVector &other) noexcept(noexcept(std::declval< BlockVector & >().storage_.swap(other.storage_)))
swap operation
Definition: bvector.hh:537
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< IsNumber< A >::value &&!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,...)
Definition: exceptions.hh:314
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:489
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:458
Dune namespace
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Implements a scalar vector view wrapper around an existing scalar.
Traits for type conversions and type information.