1 /* =========================================================================
2    Copyright (c) 2010-2016, Institute for Microelectronics,
3                             Institute for Analysis and Scientific Computing,
4                             TU Wien.
5    Portions of this software are copyright by UChicago Argonne, LLC.
6 
7                             -----------------
8                   ViennaCL - The Vienna Computing Library
9                             -----------------
10 
11    Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
12 
13    (A list of authors and contributors can be found in the PDF manual)
14 
15    License:         MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17 
18 
19 /** \file tests/src/scalar.cpp  Tests operations for viennacl::scalar objects.
20 *   \test Tests operations for viennacl::scalar objects.
21 **/
22 
23 //
24 // *** System
25 //
26 #include <iostream>
27 #include <algorithm>
28 #include <cmath>
29 
30 //
31 // *** ViennaCL
32 //
33 #include "viennacl/scalar.hpp"
34 
35 //
36 // -------------------------------------------------------------
37 //
38 template<typename ScalarType>
diff(ScalarType & s1,viennacl::scalar<ScalarType> & s2)39 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
40 {
41    viennacl::backend::finish();
42    if (std::fabs(s1 - s2) > 0)
43       return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
44    return 0;
45 }
46 //
47 // -------------------------------------------------------------
48 //
49 template< typename NumericT, typename Epsilon >
test(Epsilon const & epsilon)50 int test(Epsilon const& epsilon)
51 {
52    int retval = EXIT_SUCCESS;
53 
54    NumericT s1 = NumericT(3.1415926);
55    NumericT s2 = NumericT(2.71763);
56    NumericT s3 = NumericT(42);
57 
58    viennacl::scalar<NumericT> vcl_s1;
59    viennacl::scalar<NumericT> vcl_s2;
60    viennacl::scalar<NumericT> vcl_s3 = 1.0;
61 
62    vcl_s1 = s1;
63    if ( std::fabs(diff(s1, vcl_s1)) > epsilon )
64    {
65       std::cout << "# Error at operation: vcl_s1 = s1;" << std::endl;
66       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
67       retval = EXIT_FAILURE;
68    }
69 
70    vcl_s2 = s2;
71    if ( std::fabs(diff(s2, vcl_s2)) > epsilon )
72    {
73       std::cout << "# Error at operation: vcl_s2 = s2;" << std::endl;
74       std::cout << "  diff: " << fabs(diff(s2, vcl_s2)) << std::endl;
75       retval = EXIT_FAILURE;
76    }
77 
78    vcl_s3 = s3;
79    if ( std::fabs(diff(s3, vcl_s3)) > epsilon )
80    {
81       std::cout << "# Error at operation: vcl_s3 = s3;" << std::endl;
82       std::cout << "  diff: " << s3 - vcl_s3 << std::endl;
83       retval = EXIT_FAILURE;
84    }
85 
86    NumericT tmp = s2;
87    s2 = s1;
88    s1 = tmp;
89    viennacl::linalg::swap(vcl_s1, vcl_s2);
90    if ( fabs(diff(s1, vcl_s1)) > epsilon )
91    {
92       std::cout << "# Error at operation: swap " << std::endl;
93       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
94       retval = EXIT_FAILURE;
95    }
96 
97    s1 += s2;
98    vcl_s1 += vcl_s2;
99    if ( fabs(diff(s1, vcl_s1)) > epsilon )
100    {
101       std::cout << "# Error at operation: += " << std::endl;
102       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
103       retval = EXIT_FAILURE;
104    }
105 
106    s1 *= s3;
107    vcl_s1 *= vcl_s3;
108 
109    if ( fabs(diff(s1, vcl_s1)) > epsilon )
110    {
111       std::cout << "# Error at operation: *= " << std::endl;
112       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
113       retval = EXIT_FAILURE;
114    }
115 
116    s1 -= s2;
117    vcl_s1 -= vcl_s2;
118    if ( fabs(diff(s1, vcl_s1)) > epsilon )
119    {
120       std::cout << "# Error at operation: -= " << std::endl;
121       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
122       retval = EXIT_FAILURE;
123    }
124 
125    s1 /= s3;
126    vcl_s1 /= vcl_s3;
127 
128    if ( fabs(diff(s1, vcl_s1)) > epsilon )
129    {
130       std::cout << "# Error at operation: /= " << std::endl;
131       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
132       retval = EXIT_FAILURE;
133    }
134 
135    s1 = vcl_s1;
136 
137    s1 = s2 + s3;
138    vcl_s1 = vcl_s2 + vcl_s3;
139    if ( fabs(diff(s1, vcl_s1)) > epsilon )
140    {
141       std::cout << "# Error at operation: s1 = s2 + s3 " << std::endl;
142       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
143       retval = EXIT_FAILURE;
144    }
145 
146    s1 += s2 + s3;
147    vcl_s1 += vcl_s2 + vcl_s3;
148    if ( fabs(diff(s1, vcl_s1)) > epsilon )
149    {
150       std::cout << "# Error at operation: s1 += s2 + s3 " << std::endl;
151       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
152       retval = EXIT_FAILURE;
153    }
154 
155    s1 -= s2 + s3;
156    vcl_s1 -= vcl_s2 + vcl_s3;
157    if ( fabs(diff(s1, vcl_s1)) > epsilon )
158    {
159       std::cout << "# Error at operation: s1 -= s2 + s3 " << std::endl;
160       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
161       retval = EXIT_FAILURE;
162    }
163 
164    s1 = s2 - s3;
165    vcl_s1 = vcl_s2 - vcl_s3;
166    if ( fabs(diff(s1, vcl_s1)) > epsilon )
167    {
168       std::cout << "# Error at operation: s1 = s2 - s3 " << std::endl;
169       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
170       retval = EXIT_FAILURE;
171    }
172 
173    s1 += s2 - s3;
174    vcl_s1 += vcl_s2 - vcl_s3;
175    if ( fabs(diff(s1, vcl_s1)) > epsilon )
176    {
177       std::cout << "# Error at operation: s1 += s2 - s3 " << std::endl;
178       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
179       retval = EXIT_FAILURE;
180    }
181 
182    s1 -= s2 - s3;
183    vcl_s1 -= vcl_s2 - vcl_s3;
184    if ( fabs(diff(s1, vcl_s1)) > epsilon )
185    {
186       std::cout << "# Error at operation: s1 -= s2 - s3 " << std::endl;
187       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
188       retval = EXIT_FAILURE;
189    }
190 
191    s1 = s2 * s3;
192    vcl_s1 = vcl_s2 * vcl_s3;
193    if ( fabs(diff(s1, vcl_s1)) > epsilon )
194    {
195       std::cout << "# Error at operation: s1 = s2 * s3 " << std::endl;
196       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
197       retval = EXIT_FAILURE;
198    }
199 
200    s1 += s2 * s3;
201    vcl_s1 += vcl_s2 * vcl_s3;
202    if ( fabs(diff(s1, vcl_s1)) > epsilon )
203    {
204       std::cout << "# Error at operation: s1 += s2 * s3 " << std::endl;
205       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
206       retval = EXIT_FAILURE;
207    }
208 
209    s1 -= s2 * s3;
210    vcl_s1 -= vcl_s2 * vcl_s3;
211    if ( fabs(diff(s1, vcl_s1)) > epsilon )
212    {
213       std::cout << "# Error at operation: s1 -= s2 * s3 " << std::endl;
214       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
215       retval = EXIT_FAILURE;
216    }
217 
218    s1 = s2 / s3;
219    vcl_s1 = vcl_s2 / vcl_s3;
220    if ( fabs(diff(s1, vcl_s1)) > epsilon )
221    {
222       std::cout << "# Error at operation: s1 = s2 / s3 " << std::endl;
223       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
224       retval = EXIT_FAILURE;
225    }
226 
227    s1 += s2 / s3;
228    vcl_s1 += vcl_s2 / vcl_s3;
229    if ( fabs(diff(s1, vcl_s1)) > epsilon )
230    {
231       std::cout << "# Error at operation: s1 += s2 / s3 " << std::endl;
232       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
233       retval = EXIT_FAILURE;
234    }
235 
236    s1 -= s2 / s3;
237    vcl_s1 -= vcl_s2 / vcl_s3;
238    if ( fabs(diff(s1, vcl_s1)) > epsilon )
239    {
240       std::cout << "# Error at operation: s1 -= s2 / s3 " << std::endl;
241       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
242       retval = EXIT_FAILURE;
243    }
244 
245    // addition with factors, =
246    vcl_s1 = s1;
247 
248    s1 = s2 * s2 + s3 * s3;
249    vcl_s1 = vcl_s2 * s2 + vcl_s3 * s3;
250    if ( fabs(diff(s1, vcl_s1)) > epsilon )
251    {
252       std::cout << "# Error at operation: s1 = s2 * s2 + s3 * s3 " << std::endl;
253       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
254       retval = EXIT_FAILURE;
255    }
256    vcl_s1 = vcl_s2 * vcl_s2 + vcl_s3 * vcl_s3;
257    if ( fabs(diff(s1, vcl_s1)) > epsilon )
258    {
259       std::cout << "# Error at operation: s1 = s2 * s2 + s3 * s3, second test " << std::endl;
260       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
261       retval = EXIT_FAILURE;
262    }
263 
264    s1 = s2 * s2 + s3 / s3;
265    vcl_s1 = vcl_s2 * s2 + vcl_s3 / s3;
266    if ( fabs(diff(s1, vcl_s1)) > epsilon )
267    {
268       std::cout << "# Error at operation: s1 = s2 * s2 + s3 / s3 " << std::endl;
269       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
270       retval = EXIT_FAILURE;
271    }
272    vcl_s1 = vcl_s2 * vcl_s2 + vcl_s3 / vcl_s3;
273    if ( fabs(diff(s1, vcl_s1)) > epsilon )
274    {
275       std::cout << "# Error at operation: s1 = s2 * s2 + s3 / s3, second test " << std::endl;
276       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
277       retval = EXIT_FAILURE;
278    }
279 
280    s1 = s2 / s2 + s3 * s3;
281    vcl_s1 = vcl_s2 / s2 + vcl_s3 * s3;
282    if ( fabs(diff(s1, vcl_s1)) > epsilon )
283    {
284       std::cout << "# Error at operation: s1 = s2 / s2 + s3 * s3 " << std::endl;
285       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
286       retval = EXIT_FAILURE;
287    }
288    vcl_s1 = vcl_s2 / vcl_s2 + vcl_s3 * vcl_s3;
289    if ( fabs(diff(s1, vcl_s1)) > epsilon )
290    {
291       std::cout << "# Error at operation: s1 = s2 / s2 + s3 * s3, second test " << std::endl;
292       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
293       retval = EXIT_FAILURE;
294    }
295 
296    s1 = s2 / s2 + s3 / s3;
297    vcl_s1 = vcl_s2 / s2 + vcl_s3 / s3;
298    if ( fabs(diff(s1, vcl_s1)) > epsilon )
299    {
300       std::cout << "# Error at operation: s1 = s2 / s2 + s3 / s3 " << std::endl;
301       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
302       retval = EXIT_FAILURE;
303    }
304    vcl_s1 = vcl_s2 / vcl_s2 + vcl_s3 / vcl_s3;
305    if ( fabs(diff(s1, vcl_s1)) > epsilon )
306    {
307       std::cout << "# Error at operation: s1 = s2 / s2 + s3 / s3, second test " << std::endl;
308       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
309       retval = EXIT_FAILURE;
310    }
311 
312    // addition with factors, +=
313    vcl_s1 = s1;
314 
315    s1 += s2 * s2 + s3 * s3;
316    vcl_s1 += vcl_s2 * s2 + vcl_s3 * s3;
317    if ( fabs(diff(s1, vcl_s1)) > epsilon )
318    {
319       std::cout << "# Error at operation: s1 += s2 * s2 + s3 * s3 " << std::endl;
320       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
321       retval = EXIT_FAILURE;
322    }
323 
324    s1 += s2 * s2 + s3 / s3;
325    vcl_s1 += vcl_s2 * s2 + vcl_s3 / s3;
326    if ( fabs(diff(s1, vcl_s1)) > epsilon )
327    {
328       std::cout << "# Error at operation: s1 += s2 * s2 + s3 / s3 " << std::endl;
329       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
330       retval = EXIT_FAILURE;
331    }
332 
333    s1 += s2 / s2 + s3 * s3;
334    vcl_s1 += vcl_s2 / s2 + vcl_s3 * s3;
335    if ( fabs(diff(s1, vcl_s1)) > epsilon )
336    {
337       std::cout << "# Error at operation: s1 += s2 / s2 + s3 * s3 " << std::endl;
338       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
339       retval = EXIT_FAILURE;
340    }
341 
342    s1 += s2 / s2 + s3 / s3;
343    vcl_s1 += vcl_s2 / s2 + vcl_s3 / s3;
344    if ( fabs(diff(s1, vcl_s1)) > epsilon )
345    {
346       std::cout << "# Error at operation: s1 += s2 / s2 + s3 / s3 " << std::endl;
347       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
348       retval = EXIT_FAILURE;
349    }
350 
351    // addition with factors, -=
352    vcl_s1 = s1;
353 
354    s1 -= s2 * s2 + s3 * s3;
355    vcl_s1 -= vcl_s2 * s2 + vcl_s3 * s3;
356    if ( fabs(diff(s1, vcl_s1)) > epsilon )
357    {
358       std::cout << "# Error at operation: s1 -= s2 * s2 + s3 * s3 " << std::endl;
359       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
360       retval = EXIT_FAILURE;
361    }
362 
363    s1 -= s2 * s2 + s3 / s3;
364    vcl_s1 -= vcl_s2 * s2 + vcl_s3 / s3;
365    if ( fabs(diff(s1, vcl_s1)) > epsilon )
366    {
367       std::cout << "# Error at operation: s1 -= s2 * s2 + s3 / s3 " << std::endl;
368       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
369       retval = EXIT_FAILURE;
370    }
371 
372    s1 -= s2 / s2 + s3 * s3;
373    vcl_s1 -= vcl_s2 / s2 + vcl_s3 * s3;
374    if ( fabs(diff(s1, vcl_s1)) > epsilon )
375    {
376       std::cout << "# Error at operation: s1 -= s2 / s2 + s3 * s3 " << std::endl;
377       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
378       retval = EXIT_FAILURE;
379    }
380 
381    s1 -= s2 / s2 + s3 / s3;
382    vcl_s1 -= vcl_s2 / s2 + vcl_s3 / s3;
383    if ( fabs(diff(s1, vcl_s1)) > epsilon )
384    {
385       std::cout << "# Error at operation: s1 -= s2 / s2 + s3 / s3 " << std::endl;
386       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
387       retval = EXIT_FAILURE;
388    }
389 
390 
391    // lenghty expression:
392 
393    s1 = s2 + s3 * s2 - s3 / s1;
394    vcl_s1 = vcl_s2 + vcl_s3 * vcl_s2 - vcl_s3 / vcl_s1;
395    if ( fabs(diff(s1, vcl_s1)) > epsilon )
396    {
397       std::cout << "# Error at operation: + * - / " << std::endl;
398       std::cout << "  diff: " << fabs(diff(s1, vcl_s1)) << std::endl;
399       retval = EXIT_FAILURE;
400    }
401 
402    return retval;
403 }
404 //
405 // -------------------------------------------------------------
406 //
main()407 int main()
408 {
409    std::cout << std::endl;
410    std::cout << "----------------------------------------------" << std::endl;
411    std::cout << "----------------------------------------------" << std::endl;
412    std::cout << "## Test :: Scalar" << std::endl;
413    std::cout << "----------------------------------------------" << std::endl;
414    std::cout << "----------------------------------------------" << std::endl;
415    std::cout << std::endl;
416 
417    int retval = EXIT_SUCCESS;
418 
419    std::cout << std::endl;
420    std::cout << "----------------------------------------------" << std::endl;
421    std::cout << std::endl;
422    {
423       typedef float NumericT;
424       NumericT epsilon = NumericT(1.0E-5);
425       std::cout << "# Testing setup:" << std::endl;
426       std::cout << "  eps:     " << epsilon << std::endl;
427       std::cout << "  numeric: float" << std::endl;
428       retval = test<NumericT>(epsilon);
429       if ( retval == EXIT_SUCCESS )
430          std::cout << "# Test passed" << std::endl;
431       else
432          return retval;
433    }
434    std::cout << std::endl;
435    std::cout << "----------------------------------------------" << std::endl;
436    std::cout << std::endl;
437 #ifdef VIENNACL_WITH_OPENCL
438    if ( viennacl::ocl::current_device().double_support() )
439 #endif
440    {
441       {
442          typedef double NumericT;
443          NumericT epsilon = 1.0E-10;
444          std::cout << "# Testing setup:" << std::endl;
445          std::cout << "  eps:     " << epsilon << std::endl;
446          std::cout << "  numeric: double" << std::endl;
447          retval = test<NumericT>(epsilon);
448          if ( retval == EXIT_SUCCESS )
449            std::cout << "# Test passed" << std::endl;
450          else
451            return retval;
452       }
453       std::cout << std::endl;
454    }
455 
456   std::cout << std::endl;
457   std::cout << "------- Test completed --------" << std::endl;
458   std::cout << std::endl;
459 
460 
461    return retval;
462 }
463 //
464 // -------------------------------------------------------------
465 //
466 
467