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