1 //
2 // Copyright Toon Knapen, Karl Meerbergen
3 //
4 // Distributed under the Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //
8
9 #include "random.hpp"
10
11 #include <boost/numeric/bindings/ublas/vector_proxy.hpp>
12 #include <boost/numeric/bindings/ublas/matrix_proxy.hpp>
13 #include <boost/numeric/ublas/blas.hpp>
14 #include <boost/numeric/bindings/ublas/vector.hpp>
15 #include <boost/numeric/bindings/ublas/matrix.hpp>
16 #include <boost/numeric/bindings/std/vector.hpp>
17 #include <boost/numeric/bindings/blas/level1/axpy.hpp>
18 #include <boost/numeric/bindings/blas/level1/asum.hpp>
19 #include <boost/numeric/bindings/blas/level1/copy.hpp>
20 #include <boost/numeric/bindings/blas/level1/dot.hpp>
21 #include <boost/numeric/bindings/blas/level1/dotc.hpp>
22 #include <boost/numeric/bindings/blas/level1/dotu.hpp>
23 #include <boost/numeric/bindings/blas/level1/nrm2.hpp>
24 #include <boost/numeric/bindings/blas/level1/scal.hpp>
25 #include <boost/type_traits/is_complex.hpp>
26 #include <boost/mpl/if.hpp>
27
28 #include <vector>
29 #include <complex>
30 #include <iostream>
31 #include <limits>
32 #include <cmath>
33
34 namespace bindings = boost::numeric::bindings;
35
36 // Randomize a vector (using functions from random.hpp)
37 template <typename V>
randomize(V & v)38 void randomize(V& v)
39 {
40 for(typename V::size_type i=0; i<v.size(); ++i)
41 v[i] = random_value< typename V::value_type >() ;
42 } // randomize()
43
44
abs_sum_value(float const & f)45 float abs_sum_value(float const& f)
46 {
47 using namespace std ;
48 return abs(f) ;
49 }
50
abs_sum_value(double const & f)51 double abs_sum_value(double const& f)
52 {
53 using namespace std ;
54 return abs(f) ;
55 }
56
abs_sum_value(std::complex<float> const & f)57 float abs_sum_value(std::complex< float > const& f)
58 {
59 using namespace std ;
60 return abs(f.real()) + abs(f.imag()) ;
61 }
62
abs_sum_value(std::complex<double> const & f)63 double abs_sum_value(std::complex< double > const& f)
64 {
65 using namespace std ;
66 return abs(f.real()) + abs(f.imag()) ;
67 }
68
69 template <typename V>
abs_sum(V const & v)70 typename bindings::remove_imaginary<typename V::value_type>::type abs_sum(V const& v)
71 {
72 typedef typename bindings::remove_imaginary<typename V::value_type>::type real_type ;
73
74 real_type sum(0.0) ;
75 for(typename V::size_type i=0; i<v.size(); ++i)
76 {
77 sum += abs_sum_value(v[i]) ;
78 }
79 return sum ;
80 }
81
82
83 // Blas operations using one vector.
84 template <typename T>
85 struct OneVector
86 {
87 boost::numeric::ublas::vector<T> v_ref_ ;
88
89 // Initialize : set reference vector (ublas)
OneVectorOneVector90 OneVector()
91 : v_ref_(10)
92 {
93 randomize(v_ref_);
94 }
95
96 template <typename V>
operator ()OneVector97 int operator()(V& v) const
98 {
99 using namespace boost::numeric::bindings::blas ;
100
101 typedef typename V::value_type value_type ;
102 typedef typename bindings::remove_imaginary<value_type>::type real_type ;
103
104 // Copy vector from reference
105 for(typename V::size_type i=0; i<v_ref_.size(); ++i)
106 v[i] = v_ref_(i);
107
108 // Test blas routines and compare with reference
109 real_type nrm = nrm2(v);
110 if(std::abs(nrm - norm_2(v_ref_)) > std::numeric_limits< real_type >::epsilon() * norm_2(v_ref_))
111 {
112 std::cout << "nrm2 : " << std::abs(nrm - norm_2(v_ref_)) << " > " << std::numeric_limits< real_type >::epsilon() * norm_2(v_ref_) << std::endl ;
113 return 255 ;
114 }
115
116 nrm = asum(v);
117 if(std::abs(nrm - abs_sum(v_ref_)) > std::numeric_limits< real_type >::epsilon() * abs_sum(v_ref_))
118 {
119 std::cout << "asum : " << std::abs(nrm - abs_sum(v_ref_)) << " > " << std::numeric_limits< real_type >::epsilon() * abs_sum(v_ref_) << std::endl ;
120 return 255 ;
121 }
122
123 scal(value_type(2.0), v);
124 for(typename V::size_type i=0; i<v_ref_.size(); ++i)
125 if(std::abs(v[i] - real_type(2.0)*v_ref_(i)) > real_type(2.0)*std::abs(v_ref_(i))) return 255 ;
126
127 return 0;
128 }
129
130 // Return the size of a vector.
sizeOneVector131 size_t size() const
132 {
133 return v_ref_.size();
134 }
135 };
136
137
138 // Operations with two vectors.
139 template <typename T, typename V>
140 struct BaseTwoVectorOperations
141 {
142 typedef T value_type ;
143 typedef typename bindings::remove_imaginary<value_type>::type real_type ;
144 typedef boost::numeric::ublas::vector<T> ref_vector_type ;
145
146 // Initialize: select the first vector and set the reference vectors (ublas)
BaseTwoVectorOperationsBaseTwoVectorOperations147 BaseTwoVectorOperations(V& v, const ref_vector_type& v1_ref, const ref_vector_type& v2_ref)
148 : v_(v)
149 , v1_ref_(v1_ref)
150 , v2_ref_(v2_ref)
151 {}
152
153 // Copy the 2nd reference vector into w.
154 template <typename W>
copy_vectorBaseTwoVectorOperations155 void copy_vector(W& w) const
156 {
157 for(size_t i=0; i<size(); ++i)
158 {
159 w[i] = v2_ref_(i);
160 }
161 } // copy_vector()
162
163 // Get the size of a vector.
sizeBaseTwoVectorOperations164 size_t size() const
165 {
166 return v_.size();
167 }
168
169 // Data members.
170 V& v_ ;
171 const ref_vector_type& v1_ref_, v2_ref_ ;
172 };
173
174
175 template <typename T, typename V>
176 struct TwoVectorOperations { } ;
177
178
179 template <typename V>
180 struct TwoVectorOperations< float, V>
181 : BaseTwoVectorOperations<float,V>
182 {
183 typedef typename V::value_type value_type ;
184 typedef typename bindings::remove_imaginary<value_type>::type real_type ;
185 typedef typename BaseTwoVectorOperations<float,V>::ref_vector_type ref_vector_type ;
186
TwoVectorOperationsTwoVectorOperations187 TwoVectorOperations(V& v, const ref_vector_type& v1_ref, const ref_vector_type& v2_ref)
188 : BaseTwoVectorOperations<float,V>(v, v1_ref, v2_ref)
189 {}
190
191 // Perform the tests of blas functions and compare with reference
192 template <typename W>
operator ()TwoVectorOperations193 int operator()(W& w) const
194 {
195 using namespace boost::numeric::bindings::blas ;
196 real_type safety_factor(1.5);
197
198 copy_vector(w);
199
200 // Test blas routines
201 value_type prod = dot(this->v_, w);
202 if(std::abs(prod - inner_prod(this->v1_ref_, this->v2_ref_))
203 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(prod)) return 255 ;
204
205 axpy(value_type(2.0), this->v_, w);
206 for(size_t i=0; i<this->size(); ++i)
207 if(std::abs(w[i] - (this->v2_ref_(i) + value_type(2.0)*this->v1_ref_(i)))
208 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(w[i])) return 255 ;
209
210 scal(value_type(0.0), w) ;
211 copy(this->v_, w) ;
212 for(size_t i=0; i<this->size(); ++i)
213 {
214 if(std::abs(w[i] - this->v_[i]) != 0.0) return 255 ;
215 }
216
217 return 0;
218 }
219 };
220
221
222 template <typename V>
223 struct TwoVectorOperations< double, V>
224 : BaseTwoVectorOperations<double,V>
225 {
226 typedef typename V::value_type value_type ;
227 typedef typename bindings::remove_imaginary<value_type>::type real_type ;
228 typedef typename BaseTwoVectorOperations<double,V>::ref_vector_type ref_vector_type ;
229
TwoVectorOperationsTwoVectorOperations230 TwoVectorOperations(V& v, const ref_vector_type& v1_ref, const ref_vector_type& v2_ref)
231 : BaseTwoVectorOperations<double,V>(v, v1_ref, v2_ref)
232 {}
233
234 // Perform the tests of blas functions and compare with reference
235 template <typename W>
operator ()TwoVectorOperations236 int operator()(W& w) const
237 {
238 using namespace boost::numeric::bindings::blas ;
239
240 copy_vector(w);
241
242 // Test blas routines
243 value_type prod = dot(this->v_, w);
244 if(std::abs(prod - inner_prod(this->v1_ref_, this->v2_ref_))
245 > std::numeric_limits< real_type >::epsilon() * std::abs(prod)) return 255 ;
246
247 axpy(value_type(2.0), this->v_, w);
248 for(size_t i=0; i<this->size(); ++i)
249 if(std::abs(w[i] - (this->v2_ref_(i) + value_type(2.0)*this->v1_ref_(i)))
250 > std::numeric_limits< real_type >::epsilon() * std::abs(w[i])) return 255 ;
251
252 copy_vector(w) ;
253 scal(value_type(-1.0), w) ;
254 ::boost::numeric::bindings::blas::copy(this->v_, w) ;
255 for(size_t i=0; i<this->size(); ++i)
256 {
257 if(w[i] != this->v_[i]) return 255 ;
258 }
259
260 return 0;
261 }
262 };
263
264
265 template <typename V>
266 struct TwoVectorOperations< std::complex<float>, V>
267 : BaseTwoVectorOperations< std::complex<float>, V>
268 {
269 typedef typename V::value_type value_type ;
270 typedef typename bindings::remove_imaginary<value_type>::type real_type ;
271 typedef typename BaseTwoVectorOperations<std::complex<float>,V>::ref_vector_type ref_vector_type ;
272
TwoVectorOperationsTwoVectorOperations273 TwoVectorOperations(V& v, const ref_vector_type& v1_ref, const ref_vector_type& v2_ref)
274 : BaseTwoVectorOperations< std::complex<float>, V>(v, v1_ref, v2_ref)
275 {}
276
277 // Perform the tests of blas functions and compare with reference
278 template <typename W>
operator ()TwoVectorOperations279 int operator()(W& w) const
280 {
281 using namespace boost::numeric::bindings::blas ;
282 real_type safety_factor(1.5);
283
284 copy_vector(w);
285
286 // Test blas routines
287 value_type prod = dotc(this->v_, w);
288 if(std::abs(prod - inner_prod(conj(this->v1_ref_), this->v2_ref_))
289 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(prod)) return 255 ;
290
291 prod = dotu(this->v_, w);
292 if(std::abs(prod - inner_prod(this->v1_ref_, this->v2_ref_))
293 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(prod)) return 255 ;
294
295 axpy(value_type(2.0), this->v_, w);
296 for(size_t i=0; i<this->size(); ++i)
297 if(std::abs(w[i] - (this->v2_ref_(i) + value_type(2.0)*this->v1_ref_(i)))
298 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(w[i])) return 255 ;
299
300 scal(value_type(0.0), w) ;
301 copy(this->v_, w) ;
302 for(size_t i=0; i<this->size(); ++i)
303 {
304 if(std::abs(w[i] - this->v_[i]) != 0.0) return 255 ;
305 }
306
307 return 0;
308 }
309 };
310
311
312 template <typename V>
313 struct TwoVectorOperations< std::complex<double>, V>
314 : BaseTwoVectorOperations< std::complex<double>, V>
315 {
316 typedef typename V::value_type value_type ;
317 typedef typename bindings::remove_imaginary<value_type>::type real_type ;
318 typedef typename BaseTwoVectorOperations<std::complex<double>,V>::ref_vector_type ref_vector_type ;
319
TwoVectorOperationsTwoVectorOperations320 TwoVectorOperations(V& v, const ref_vector_type& v1_ref, const ref_vector_type& v2_ref)
321 : BaseTwoVectorOperations< std::complex<double>, V>(v, v1_ref, v2_ref)
322 {}
323
324 // Perform the tests of blas functions and compare with reference
325 template <typename W>
operator ()TwoVectorOperations326 int operator()(W& w) const
327 {
328 using namespace boost::numeric::bindings::blas ;
329 real_type safety_factor(1.5);
330
331 copy_vector(w);
332
333 // Test blas routines
334 value_type prod = dotc(this->v_, w);
335 if(std::abs(prod - inner_prod(conj(this->v1_ref_), this->v2_ref_))
336 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(prod)) return 255 ;
337
338 prod = dotu(this->v_, w);
339 if(std::abs(prod - inner_prod(this->v1_ref_, this->v2_ref_))
340 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(prod)) return 255 ;
341
342 axpy(value_type(2.0), this->v_, w);
343 for(size_t i=0; i<this->size(); ++i)
344 if(std::abs(w[i] - (this->v2_ref_(i) + value_type(2.0)*this->v1_ref_(i)))
345 > safety_factor*std::numeric_limits< real_type >::epsilon() * std::abs(w[i])) return 255 ;
346
347 scal(value_type(0.0), w) ;
348 copy(this->v_, w) ;
349 for(size_t i=0; i<this->size(); ++i)
350 {
351 if(std::abs(w[i] - this->v_[i]) != 0.0) return 255 ;
352 }
353
354 return 0;
355 }
356 };
357
358
359 // Run the tests for different types of vectors.
360 template <typename T, typename F>
different_vectors(const F & f)361 int different_vectors(const F& f)
362 {
363 // Do test for different types of vectors
364 {
365 std::cout << " ublas::vector\n" ;
366 boost::numeric::ublas::vector< T > v(f.size());
367 if(f(v)) return 255 ;
368 }
369 {
370 std::cout << " std::vector\n" ;
371 std::vector<T> v_ref(f.size());
372 if(f(v_ref)) return 255 ;
373 }
374 {
375 std::cout << " ublas::vector_range\n" ;
376 typedef boost::numeric::ublas::vector< T > vector_type ;
377 vector_type v(f.size()*2);
378 boost::numeric::ublas::vector_range< vector_type > vr(v, boost::numeric::ublas::range(1,1+f.size()));
379 if(f(vr)) return 255 ;
380 }
381 {
382 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > matrix_type ;
383 matrix_type m(f.size(),f.size()) ;
384
385 std::cout << " ublas::matrix_column\n" ;
386 boost::numeric::ublas::matrix_column< matrix_type > m_c(m, 2);
387 if(f(m_c)) return 255 ;
388
389 std::cout << " ublas::matrix_row\n" ;
390 boost::numeric::ublas::matrix_row< matrix_type > m_r(m, 1);
391 if(f(m_r)) return 255 ;
392 }
393 return 0;
394 } // different_vectors()
395
396
397 // This is the functor that selects the first vector of the tests that use two vectors.
398 template <typename T>
399 struct TwoVector
400 {
TwoVectorTwoVector401 TwoVector()
402 : v1_ref_(10)
403 , v2_ref_(10)
404 {}
405
406 template <typename V>
operator ()TwoVector407 int operator()(V& v) const
408 {
409 for(size_t i=0; i<size(); ++i) v[i] = v1_ref_(i) ;
410 return different_vectors<T,TwoVectorOperations<T,V> >(TwoVectorOperations<T,V>(v, v1_ref_, v2_ref_)) ;
411 }
412
sizeTwoVector413 size_t size() const
414 {
415 return v1_ref_.size() ;
416 }
417
418 boost::numeric::ublas::vector<T> v1_ref_ ;
419 boost::numeric::ublas::vector<T> v2_ref_ ;
420 }; // TwoVector
421
422
423 // Run the test for a specific value_type T.
424 template <typename T>
do_value_type()425 int do_value_type()
426 {
427 // Tests for functions with one vector argument.
428 std::cout << " one argument\n";
429 if(different_vectors<T,OneVector<T> >(OneVector<T> ())) return 255 ;
430
431 // Tests for functions with two vector arguments.
432 std::cout << " two arguments\n";
433 if(different_vectors<T,TwoVector<T> >(TwoVector<T>())) return 255;
434 return 0;
435 } // do_value_type()
436
437
main()438 int main()
439 {
440 // Run regression for Real/Complex
441 std::cout << "float\n";
442 if(do_value_type<float>()) return 255 ;
443 std::cout << "double\n";
444 if(do_value_type<double>()) return 255 ;
445 std::cout << "complex<float>\n";
446 if(do_value_type<std::complex<float> >()) return 255 ;
447 std::cout << "complex<double>\n";
448 if(do_value_type<std::complex<double> >()) return 255 ;
449
450 std::cout << "Regression test successful\n" ;
451
452 return 0 ;
453 }
454
455
456