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])18void 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()72bool 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 * [])95int 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