123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103 |
- /* Boost test/det.cpp
- * test protected and unprotected rounding on an unstable determinant
- *
- * Copyright 2002-2003 Guillaume Melquiond
- *
- * Distributed under the Boost Software License, Version 1.0.
- * (See accompanying file LICENSE_1_0.txt or
- * copy at http://www.boost.org/LICENSE_1_0.txt)
- */
- #include <boost/numeric/interval.hpp>
- #include <boost/test/minimal.hpp>
- #include "bugs.hpp"
- #define size 8
- template<class I>
- void det(I (&mat)[size][size]) {
- for(int i = 0; i < size; i++)
- for(int j = 0; j < size; j++)
- mat[i][j] = I(1) / I(i + j + 1);
- for(int i = 0; i < size - 1; i++) {
- int m = i, n = i;
- typename I::base_type v = 0;
- for(int a = i; a < size; a++)
- for(int b = i; b < size; b++) {
- typename I::base_type w = abs(mat[a][b]).lower();
- if (w > v) { m = a; n = b; v = w; }
- }
- if (n != i)
- for(int a = 0; a < size; a++) {
- I t = mat[a][n];
- mat[a][n] = mat[a][i];
- mat[a][i] = t;
- }
- if (m != i)
- for(int b = i; b < size; b++) {
- I t = mat[m][b];
- mat[m][b] = mat[m][i];
- mat[m][i] = t;
- }
- if (((m + n) & 1) == 1) { };
- I c = mat[i][i];
- for(int j = i + 1; j < size; j++) {
- I f = mat[j][i] / c;
- for(int k = i; k < size; k++)
- mat[j][k] -= f * mat[i][k];
- }
- if (zero_in(c)) return;
- }
- }
- namespace my_namespace {
- using namespace boost;
- using namespace numeric;
- using namespace interval_lib;
- template<class T>
- struct variants {
- typedef interval<T> I_op;
- typedef typename change_rounding<I_op, save_state<rounded_arith_std<T> > >::type I_sp;
- typedef typename unprotect<I_op>::type I_ou;
- typedef typename unprotect<I_sp>::type I_su;
- typedef T type;
- };
- }
- template<class T>
- bool test() {
- typedef my_namespace::variants<double> types;
- types::I_op mat_op[size][size];
- types::I_sp mat_sp[size][size];
- types::I_ou mat_ou[size][size];
- types::I_su mat_su[size][size];
- det(mat_op);
- det(mat_sp);
- { types::I_op::traits_type::rounding rnd; det(mat_ou); }
- { types::I_sp::traits_type::rounding rnd; det(mat_su); }
- for(int i = 0; i < size; i++)
- for(int j = 0; j < size; j++) {
- typedef types::I_op I;
- I d_op = mat_op[i][j];
- I d_sp = mat_sp[i][j];
- I d_ou = mat_ou[i][j];
- I d_su = mat_su[i][j];
- if (!(equal(d_op, d_sp) && equal(d_sp, d_ou) && equal(d_ou, d_su)))
- return false;
- }
- return true;
- }
- int test_main(int, char *[]) {
- BOOST_CHECK(test<float>());
- BOOST_CHECK(test<double>());
- BOOST_CHECK(test<long double>());
- # ifdef __BORLANDC__
- ::detail::ignore_warnings();
- # endif
- return 0;
- }
|