1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #include "config.h"
4 #include <dune/istl/bvector.hh>
5 #include <dune/istl/vbvector.hh>
6 #include <dune/common/fvector.hh>
7 #include <dune/common/classname.hh>
8 #include <dune/common/std/type_traits.hh>
9 #include <dune/common/test/testsuite.hh>
10 #include <dune/common/scalarvectorview.hh>
11
12 template<typename T>
13 struct Sign
14 {
complexSignSign15 static T complexSign()
16 { return T(1.0); }
sqrtComplexSignSign17 static T sqrtComplexSign()
18 { return T(1.0); }
19 };
20
21 template<typename T>
22 struct Sign< std::complex<T> >
23 {
complexSignSign24 static std::complex<T> complexSign()
25 { return std::complex<T>(-1.0); }
sqrtComplexSignSign26 static std::complex<T> sqrtComplexSign()
27 { return std::complex<T>(0.0, 1.0); }
28 };
29
30 // scalar ordering doesn't work for complex numbers
31 template <class RealBlockVector, class ComplexBlockVector>
DotProductTest(const size_t numBlocks,const size_t blockSizeOrCapacity)32 Dune::TestSuite DotProductTest(const size_t numBlocks,const size_t blockSizeOrCapacity) {
33 Dune::TestSuite t;
34
35 typedef typename RealBlockVector::field_type rt;
36 typedef typename ComplexBlockVector::field_type ct;
37 const rt myEps((rt)1e-6);
38
39 static_assert(std::is_same< typename Dune::FieldTraits<rt>::real_type, rt>::value,
40 "DotProductTest requires real data type for first block vector!");
41
42 const ct complexSign = Sign<ct>::complexSign();
43 const ct I = Sign<ct>::sqrtComplexSign();
44
45 typedef typename RealBlockVector::size_type size_type;
46
47 // empty vectors
48 RealBlockVector one(numBlocks,blockSizeOrCapacity);
49 ComplexBlockVector iVec(numBlocks,blockSizeOrCapacity);
50
51 using RealBlockType = typename RealBlockVector::block_type;
52
53 const size_type blockSize = Dune::Impl::asVector(one[0]).size();
54
55 t.require(numBlocks==one.N());
56 t.require(numBlocks==iVec.N());
57 const size_type length = numBlocks * blockSize; // requires inner block size of VariableBlockVector to be 1!
58
59 ct ctlength = ct(length);
60
61 std::cout << __func__ << "\t \t ( " << Dune::className(one) << " and \n \t \t \t " << Dune::className(iVec) << " )" << std::endl << std::endl;
62
63 // initialize vectors with data
64 for(size_type i=0; i < numBlocks; ++i) {
65 one[i] = 1.;
66 iVec[i] = I;
67 }
68
69 ct result = ct();
70
71 // blockwise dot tests
72 if constexpr (!Dune::IsNumber<RealBlockType>{}){
73 result = ct();
74 for(size_type i=0; i < numBlocks; ++i) {
75 result += dot(one[i],one[i]) + (one[i]).dot(one[i]);
76 }
77
78 t.check(std::abs(result-ct(2)*ctlength)<= myEps);
79
80 result = ct();
81 for(size_type i=0; i < numBlocks; ++i) {
82 result += dot(iVec[i],iVec[i])+ (iVec[i]).dot(iVec[i]);
83 }
84
85 t.check(std::abs(result-ct(2)*ctlength)<= myEps);
86
87 // blockwise dotT / operator * tests
88 result = ct();
89 for(size_type i=0; i < numBlocks; ++i) {
90 result += dotT(one[i],one[i]) + one[i]*one[i];
91 }
92
93 t.check(std::abs(result-ct(2)*ctlength)<= myEps);
94
95 result = ct();
96 for(size_type i=0; i < numBlocks; ++i) {
97 result += dotT(iVec[i],iVec[i]) + iVec[i]*iVec[i];
98 }
99
100 t.check(std::abs(result-complexSign*ct(2)*ctlength)<= myEps);
101 }
102
103 // global operator * tests
104 result = one*one + dotT(one,one);
105
106 t.check(std::abs(result-ct(2)*ctlength)<= myEps);
107
108 result = iVec*iVec + dotT(iVec,iVec);
109
110 t.check(std::abs(result-complexSign*ct(2)*ctlength)<= myEps);
111
112 // global operator dot(,) tests
113 result = one.dot(one) + dot(one,one);
114
115 t.check(std::abs(result-ct(2)*ctlength)<= myEps);
116 result = iVec.dot(iVec) + dot(iVec,iVec);
117
118 t.check(std::abs(result-ct(2)*ctlength)<= myEps);
119
120 // mixed global dotT tests
121 result = iVec*one + one*iVec + dotT(one,iVec) + dotT(iVec,one);
122
123 t.check(std::abs(result-ct(4)*ctlength*I)<= myEps);
124
125 // mixed global dot tests
126 result = iVec.dot(one) + dot(iVec,one);
127
128 t.check(std::abs(result-ct(2)*complexSign*ctlength*I)<= myEps);
129 result = one.dot(iVec) + dot(one,iVec);
130
131 t.check(std::abs(result-ct(2)*ctlength*I)<= myEps);
132
133 return t;
134 }
135
136
main()137 int main()
138 {
139 Dune::TestSuite t;
140 const size_t BlockSize = 5;
141 const size_t numBlocks = 10;
142 const size_t capacity = BlockSize * numBlocks * 2; // use capacity here, that we can use the a constructor taking two integers for both BlockVector and VariableBlockVector
143
144 t.subTest(DotProductTest<Dune::BlockVector<Dune::FieldVector<int,BlockSize> >, Dune::BlockVector<Dune::FieldVector<int,BlockSize> > > (numBlocks,capacity));
145 t.subTest(DotProductTest<Dune::VariableBlockVector<Dune::FieldVector<int,1> >, Dune::VariableBlockVector<Dune::FieldVector<int,1> > > (numBlocks,1));
146
147 t.subTest(DotProductTest<Dune::BlockVector<Dune::FieldVector<float,BlockSize> >, Dune::BlockVector<Dune::FieldVector<std::complex<float>,BlockSize> > > (numBlocks,capacity));
148 t.subTest(DotProductTest<Dune::VariableBlockVector<Dune::FieldVector<float,1> >, Dune::VariableBlockVector<Dune::FieldVector<std::complex<float>,1> > > (numBlocks,BlockSize));
149
150 t.subTest(DotProductTest<Dune::BlockVector<Dune::FieldVector<float,BlockSize> >, Dune::BlockVector<Dune::FieldVector<float,BlockSize> > > (numBlocks,capacity));
151 t.subTest(DotProductTest<Dune::VariableBlockVector<Dune::FieldVector<float,1> >, Dune::VariableBlockVector<Dune::FieldVector<float,1> > > (numBlocks,1));
152
153 t.subTest(DotProductTest<Dune::BlockVector<Dune::FieldVector<double,BlockSize> >, Dune::BlockVector<Dune::FieldVector<std::complex<double>,BlockSize> > > (numBlocks,capacity));
154 t.subTest(DotProductTest<Dune::VariableBlockVector<Dune::FieldVector<double,1> >, Dune::VariableBlockVector<Dune::FieldVector<std::complex<double>,1> > > (numBlocks,BlockSize));
155
156 t.subTest(DotProductTest<Dune::BlockVector<Dune::FieldVector<double,BlockSize> >, Dune::BlockVector<Dune::FieldVector<double,BlockSize> > > (numBlocks,capacity));
157 t.subTest(DotProductTest<Dune::VariableBlockVector<Dune::FieldVector<double,1> >, Dune::VariableBlockVector<Dune::FieldVector<double,1> > > (numBlocks,BlockSize));
158
159 t.subTest(DotProductTest<Dune::BlockVector<double>, Dune::BlockVector<double> > (numBlocks,capacity));
160
161 return t.exit();
162 }
163