| /* Boost example/filter.cpp |
| * two examples of filters for computing the sign of a determinant |
| * the second filter is based on an idea presented in |
| * "Interval arithmetic yields efficient dynamic filters for computational |
| * geometry" by Brönnimann, Burnikel and Pion, 2001 |
| * |
| * Copyright 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 <iostream> |
| |
| namespace dummy { |
| using namespace boost; |
| using namespace numeric; |
| using namespace interval_lib; |
| typedef save_state<rounded_arith_opp<double> > R; |
| typedef checking_no_nan<double, checking_no_empty<double> > P; |
| typedef interval<double, policies<R, P> > I; |
| } |
| |
| template<class T> |
| class vector { |
| T* ptr; |
| public: |
| vector(int d) { ptr = (T*)malloc(sizeof(T) * d); } |
| ~vector() { free(ptr); } |
| const T& operator[](int i) const { return ptr[i]; } |
| T& operator[](int i) { return ptr[i]; } |
| }; |
| |
| template<class T> |
| class matrix { |
| int dim; |
| T* ptr; |
| public: |
| matrix(int d): dim(d) { ptr = (T*)malloc(sizeof(T) * dim * dim); } |
| ~matrix() { free(ptr); } |
| int get_dim() const { return dim; } |
| void assign(const matrix<T> &a) { memcpy(ptr, a.ptr, sizeof(T) * dim * dim); } |
| const T* operator[](int i) const { return &(ptr[i * dim]); } |
| T* operator[](int i) { return &(ptr[i * dim]); } |
| }; |
| |
| typedef dummy::I I_dbl; |
| |
| /* compute the sign of a determinant using an interval LU-decomposition; the |
| function answers 1 or -1 if the determinant is positive or negative (and |
| more importantly, the result must be provable), or 0 if the algorithm was |
| unable to get a correct sign */ |
| int det_sign_algo1(const matrix<double> &a) { |
| int dim = a.get_dim(); |
| vector<int> p(dim); |
| for(int i = 0; i < dim; i++) p[i] = i; |
| int sig = 1; |
| I_dbl::traits_type::rounding rnd; |
| typedef boost::numeric::interval_lib::unprotect<I_dbl>::type I; |
| matrix<I> u(dim); |
| for(int i = 0; i < dim; i++) { |
| const double* line1 = a[i]; |
| I* line2 = u[i]; |
| for(int j = 0; j < dim; j++) |
| line2[j] = line1[j]; |
| } |
| // computation of L and U |
| for(int i = 0; i < dim; i++) { |
| // partial pivoting |
| { |
| int pivot = i; |
| double max = 0; |
| for(int j = i; j < dim; j++) { |
| const I &v = u[p[j]][i]; |
| if (zero_in(v)) continue; |
| double m = norm(v); |
| if (m > max) { max = m; pivot = j; } |
| } |
| if (max == 0) return 0; |
| if (pivot != i) { |
| sig = -sig; |
| int tmp = p[i]; |
| p[i] = p[pivot]; |
| p[pivot] = tmp; |
| } |
| } |
| // U[i,?] |
| { |
| I *line1 = u[p[i]]; |
| const I &pivot = line1[i]; |
| if (boost::numeric::interval_lib::cerlt(pivot, 0.)) sig = -sig; |
| for(int k = i + 1; k < dim; k++) { |
| I *line2 = u[p[k]]; |
| I fact = line2[i] / pivot; |
| for(int j = i + 1; j < dim; j++) line2[j] -= fact * line1[j]; |
| } |
| } |
| } |
| return sig; |
| } |
| |
| /* compute the sign of a determinant using a floating-point LU-decomposition |
| and an a posteriori interval validation; the meaning of the answer is the |
| same as previously */ |
| int det_sign_algo2(const matrix<double> &a) { |
| int dim = a.get_dim(); |
| vector<int> p(dim); |
| for(int i = 0; i < dim; i++) p[i] = i; |
| int sig = 1; |
| matrix<double> lui(dim); |
| { |
| // computation of L and U |
| matrix<double> lu(dim); |
| lu.assign(a); |
| for(int i = 0; i < dim; i++) { |
| // partial pivoting |
| { |
| int pivot = i; |
| double max = std::abs(lu[p[i]][i]); |
| for(int j = i + 1; j < dim; j++) { |
| double m = std::abs(lu[p[j]][i]); |
| if (m > max) { max = m; pivot = j; } |
| } |
| if (max == 0) return 0; |
| if (pivot != i) { |
| sig = -sig; |
| int tmp = p[i]; |
| p[i] = p[pivot]; |
| p[pivot] = tmp; |
| } |
| } |
| // L[?,i] and U[i,?] |
| { |
| double *line1 = lu[p[i]]; |
| double pivot = line1[i]; |
| if (pivot < 0) sig = -sig; |
| for(int k = i + 1; k < dim; k++) { |
| double *line2 = lu[p[k]]; |
| double fact = line2[i] / pivot; |
| line2[i] = fact; |
| for(int j = i + 1; j < dim; j++) line2[j] -= line1[j] * fact; |
| } |
| } |
| } |
| |
| // computation of approximate inverses: Li and Ui |
| for(int j = 0; j < dim; j++) { |
| for(int i = j + 1; i < dim; i++) { |
| double *line = lu[p[i]]; |
| double s = - line[j]; |
| for(int k = j + 1; k < i; k++) s -= line[k] * lui[k][j]; |
| lui[i][j] = s; |
| } |
| lui[j][j] = 1 / lu[p[j]][j]; |
| for(int i = j - 1; i >= 0; i--) { |
| double *line = lu[p[i]]; |
| double s = 0; |
| for(int k = i + 1; k <= j; k++) s -= line[k] * lui[k][j]; |
| lui[i][j] = s / line[i]; |
| } |
| } |
| } |
| |
| // norm of PAUiLi-I computed with intervals |
| { |
| I_dbl::traits_type::rounding rnd; |
| typedef boost::numeric::interval_lib::unprotect<I_dbl>::type I; |
| vector<I> m1(dim); |
| vector<I> m2(dim); |
| for(int i = 0; i < dim; i++) { |
| for(int j = 0; j < dim; j++) m1[j] = 0; |
| const double *l1 = a[p[i]]; |
| for(int j = 0; j < dim; j++) { |
| double v = l1[j]; // PA[i,j] |
| double *l2 = lui[j]; // Ui[j,?] |
| for(int k = j; k < dim; k++) { |
| using boost::numeric::interval_lib::mul; |
| m1[k] += mul<I>(v, l2[k]); // PAUi[i,k] |
| } |
| } |
| for(int j = 0; j < dim; j++) m2[j] = m1[j]; // PAUi[i,j] * Li[j,j] |
| for(int j = 1; j < dim; j++) { |
| const I &v = m1[j]; // PAUi[i,j] |
| double *l2 = lui[j]; // Li[j,?] |
| for(int k = 0; k < j; k++) |
| m2[k] += v * l2[k]; // PAUiLi[i,k] |
| } |
| m2[i] -= 1; // PAUiLi-I |
| double ss = 0; |
| for(int i = 0; i < dim; i++) ss = rnd.add_up(ss, norm(m2[i])); |
| if (ss >= 1) return 0; |
| } |
| } |
| return sig; |
| } |
| |
| int main() { |
| int dim = 20; |
| matrix<double> m(dim); |
| for(int i = 0; i < dim; i++) for(int j = 0; j < dim; j++) |
| m[i][j] = /*1 / (i-j-0.001)*/ cos(1+i*sin(1 + j)) /*1./(1+i+j)*/; |
| |
| // compute the sign of the determinant of a "strange" matrix with the two |
| // algorithms, the first should fail and the second succeed |
| std::cout << det_sign_algo1(m) << " " << det_sign_algo2(m) << std::endl; |
| } |