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
5 // hack to ensure assert() does something
6 // really, assert() should not be used in unit tests.
7 #ifdef NDEBUG
8 #undef NDEBUG
9 #endif
10
11 #include <cassert>
12 #include <cmath>
13 #include <complex>
14 #include <cstdlib>
15 #include <iostream>
16 #include <memory>
17
18 #include <dune/common/classname.hh>
19 #if HAVE_MPROTECT
20 #include <dune/common/debugallocator.hh>
21 #endif
22 #include <dune/common/fvector.hh>
23 #include <dune/common/poolallocator.hh>
24 #include <dune/common/scalarvectorview.hh>
25
26 #include <dune/istl/bvector.hh>
27 #include <dune/istl/test/vectortest.hh>
28
29 template<typename VectorBlock, typename V>
assign(VectorBlock & b,const V & i)30 void assign(VectorBlock& b, const V& i)
31 {
32 for (auto& entry : Dune::Impl::asVector(b))
33 entry = i;
34 }
35
36
37 template<class VectorBlock, class A=std::allocator<void> >
testVector()38 int testVector()
39 {
40 using Alloc = typename std::allocator_traits<A>::template rebind_alloc<VectorBlock>;
41 typedef Dune::BlockVector<VectorBlock, Alloc> Vector;
42
43 // empty vector
44 Vector v, w, v1(20), v2(20,100);
45 #ifdef FAIL
46 Vector v3(20,100.0);
47 #endif
48 v.reserve(100);
49 assert(100==v.capacity());
50 assert(20==v1.capacity());
51 assert(100==v2.capacity());
52 assert(20==v1.N());
53 assert(20==v2.N());
54
55 v.resize(25);
56
57 assert(25==v.N());
58
59 for(typename Vector::size_type i=0; i < v.N(); ++i)
60 v[i] = i;
61
62 for(typename Vector::size_type i=0; i < v2.N(); ++i)
63 v2[i] = i*10;
64 w = v;
65
66 testHomogeneousRandomAccessContainer(v);
67 Dune::testConstructibility<Vector>();
68 testNorms(v);
69 testVectorSpaceOperations(v);
70 testScalarProduct(v);
71
72 assert(w.N()==v.N());
73
74 for(typename Vector::size_type i=0; i < v.N(); ++i)
75 assert(v[i] == w[i]);
76
77 Vector z(w);
78
79 assert(w.N()==z.N());
80 assert(w.capacity()==z.capacity());
81
82 for(typename Vector::size_type i=0; i < w.N(); ++i)
83 assert(z[i] == w[i]);
84
85 v.reserve(150);
86 assert(150==v.capacity());
87 assert(25==v.N());
88
89 VectorBlock b;
90
91 // check the entries
92 for(typename Vector::size_type i=0; i < v.N(); ++i) {
93 assign(b, i);
94 assert(v[i] == b);
95 }
96
97 return 0;
98 }
99
testCapacity()100 void testCapacity()
101 {
102 typedef Dune::FieldVector<double,2> SmallVector;
103 typedef Dune::BlockVector<Dune::BlockVector<SmallVector> > ThreeLevelVector;
104 ThreeLevelVector vec;
105 vec.reserve(10);
106 vec.resize(10);
107 for(int i=0; i<10; ++i)
108 vec[i]=Dune::BlockVector<SmallVector>(10);
109 ThreeLevelVector vec1=vec;
110 }
111
112 template <class V>
checkNormNAN(V const & v,int line)113 void checkNormNAN(V const &v, int line) {
114 if (!std::isnan(v.one_norm())) {
115 std::cerr << "error: norm not NaN: one_norm() on line "
116 << line << " (type: " << Dune::className(v[0][0]) << ")"
117 << std::endl;
118 std::exit(-1);
119 }
120 if (!std::isnan(v.two_norm())) {
121 std::cerr << "error: norm not NaN: two_norm() on line "
122 << line << " (type: " << Dune::className(v[0][0]) << ")"
123 << std::endl;
124 std::exit(-1);
125 }
126 if (!std::isnan(v.infinity_norm())) {
127 std::cerr << "error: norm not NaN: infinity_norm() on line "
128 << line << " (type: " << Dune::className(v[0][0]) << ")"
129 << std::endl;
130 std::exit(-1);
131 }
132 }
133
134 // Make sure that vectors with NaN entries have norm NaN.
135 // See also bug flyspray/FS#1147
136 template <typename T>
137 void
test_nan(T const & mynan)138 test_nan(T const &mynan)
139 {
140 using FV = Dune::FieldVector<T,2>;
141 using V = Dune::BlockVector<FV>;
142 T n(0);
143 {
144 V v = {
145 { mynan, n },
146 { n, n }
147 };
148 checkNormNAN(v, __LINE__);
149 }
150 {
151 V v = {
152 { n, mynan },
153 { n, n }
154 };
155 checkNormNAN(v, __LINE__);
156 }
157 {
158 V v = {
159 { n, n },
160 { mynan, n }
161 };
162 checkNormNAN(v, __LINE__);
163 }
164 {
165 V v = {
166 { n, n },
167 { n, mynan }
168 };
169 checkNormNAN(v, __LINE__);
170 }
171 {
172 V v = {
173 { mynan, mynan },
174 { mynan, mynan }
175 };
176 checkNormNAN(v, __LINE__);
177 }
178 }
179
main()180 int main()
181 {
182 typedef std::complex<double> value_type;
183 //typedef double value_type;
184 typedef Dune::FieldVector<value_type,1> VectorBlock;
185 typedef Dune::BlockVector<VectorBlock> Vector;
186 Vector v;
187 v=0;
188 Dune::BlockVector<Dune::FieldVector<std::complex<double>,1> > v1;
189 v1=0;
190
191 // Test a BlockVector of BlockVectors
192 typedef Dune::BlockVector<Vector> VectorOfVector;
193 VectorOfVector vv = {{1.0, 2.0}, {3.0, 4.0, 5.0}, {6.0}};
194
195 testHomogeneousRandomAccessContainer(vv);
196 Dune::testConstructibility<VectorOfVector>();
197 testNorms(vv);
198 testVectorSpaceOperations(vv);
199 testScalarProduct(vv);
200
201 // Test construction from initializer_list
202 Vector fromInitializerList = {0,1,2};
203 assert(fromInitializerList.size() == 3);
204 assert(fromInitializerList[0] == value_type(0));
205 assert(fromInitializerList[1] == value_type(1));
206 assert(fromInitializerList[2] == value_type(2));
207
208 {
209 double nan = std::nan("");
210 test_nan(nan);
211 }
212 {
213 std::complex<double> nan( std::nan(""), 17 );
214 test_nan(nan);
215 }
216
217 int ret = 0;
218
219 ret += testVector<Dune::FieldVector<double,1> >();
220 // ret += testVector<1, Dune::PoolAllocator<void,1000000> >();
221 #if HAVE_MPROTECT
222 ret += testVector<Dune::FieldVector<double,1> , Dune::DebugAllocator<void> >();
223 #endif
224 ret += testVector<Dune::FieldVector<double,3> >();
225 // ret += testVector<3, Dune::PoolAllocator<void,1000000> >();
226 #if HAVE_MPROTECT
227 ret += testVector<Dune::FieldVector<double,1> , Dune::DebugAllocator<void> >();
228 #endif
229
230 ret += testVector<double>();
231 ret += testVector<std::complex<double> >();
232
233 testCapacity();
234
235 return ret;
236 }
237