1 /* Boost test/det.cpp
2  * test protected and unprotected rounding on an unstable determinant
3  *
4  * Copyright 2002-2003 Guillaume Melquiond
5  *
6  * Distributed under the Boost Software License, Version 1.0.
7  * (See accompanying file LICENSE_1_0.txt or
8  * copy at http://www.boost.org/LICENSE_1_0.txt)
9  */
10 
11 #include <boost/numeric/interval.hpp>
12 #include <boost/test/minimal.hpp>
13 #include "bugs.hpp"
14 
15 #define size 8
16 
17 template<class I>
det(I (& mat)[size][size])18 void det(I (&mat)[size][size]) {
19   for(int i = 0; i < size; i++)
20     for(int j = 0; j < size; j++)
21       mat[i][j] = I(1) / I(i + j + 1);
22 
23   for(int i = 0; i < size - 1; i++) {
24     int m = i, n = i;
25     typename I::base_type v = 0;
26     for(int a = i; a < size; a++)
27       for(int b = i; b < size; b++) {
28         typename I::base_type  w = abs(mat[a][b]).lower();
29         if (w > v) { m = a; n = b; v = w; }
30       }
31     if (n != i)
32       for(int a = 0; a < size; a++) {
33         I t = mat[a][n];
34         mat[a][n] = mat[a][i];
35         mat[a][i] = t;
36       }
37     if (m != i)
38       for(int b = i; b < size; b++) {
39         I t = mat[m][b];
40         mat[m][b] = mat[m][i];
41         mat[m][i] = t;
42       }
43     if (((m + n) & 1) == 1) { };
44     I c = mat[i][i];
45     for(int j = i + 1; j < size; j++) {
46       I f = mat[j][i] / c;
47       for(int k = i; k < size; k++)
48         mat[j][k] -= f * mat[i][k];
49     }
50     if (zero_in(c)) return;
51   }
52 }
53 
54 namespace my_namespace {
55 
56 using namespace boost;
57 using namespace numeric;
58 using namespace interval_lib;
59 
60 template<class T>
61 struct variants {
62   typedef interval<T> I_op;
63   typedef typename change_rounding<I_op, save_state<rounded_arith_std<T> > >::type I_sp;
64   typedef typename unprotect<I_op>::type I_ou;
65   typedef typename unprotect<I_sp>::type I_su;
66   typedef T type;
67 };
68 
69 }
70 
71 template<class T>
test()72 bool test() {
73   typedef my_namespace::variants<double> types;
74   types::I_op mat_op[size][size];
75   types::I_sp mat_sp[size][size];
76   types::I_ou mat_ou[size][size];
77   types::I_su mat_su[size][size];
78   det(mat_op);
79   det(mat_sp);
80   { types::I_op::traits_type::rounding rnd; det(mat_ou); }
81   { types::I_sp::traits_type::rounding rnd; det(mat_su); }
82   for(int i = 0; i < size; i++)
83     for(int j = 0; j < size; j++) {
84       typedef types::I_op I;
85       I d_op = mat_op[i][j];
86       I d_sp = mat_sp[i][j];
87       I d_ou = mat_ou[i][j];
88       I d_su = mat_su[i][j];
89       if (!(equal(d_op, d_sp) && equal(d_sp, d_ou) && equal(d_ou, d_su)))
90         return false;
91     }
92   return true;
93 }
94 
test_main(int,char * [])95 int test_main(int, char *[]) {
96   BOOST_CHECK(test<float>());
97   BOOST_CHECK(test<double>());
98   BOOST_CHECK(test<long double>());
99 # ifdef __BORLANDC__
100   ::detail::ignore_warnings();
101 # endif
102   return 0;
103 }
104