ceph/src/boost/libs/numeric/interval/examples/filter.cpp

   1 /* Boost example/filter.cpp
   2  * two examples of filters for computing the sign of a determinant
   3  * the second filter is based on an idea presented in
   4  * "Interval arithmetic yields efficient dynamic filters for computational
   5  * geometry" by Brönnimann, Burnikel and Pion, 2001
   6  *
   7  * Copyright 2003 Guillaume Melquiond
   8  *
   9  * Distributed under the Boost Software License, Version 1.0.
  10  * (See accompanying file LICENSE_1_0.txt or
  11  * copy at http://www.boost.org/LICENSE_1_0.txt)
  12  */
  13
  14 #include <boost/numeric/interval.hpp>
  15 #include <iostream>
  16
  17 namespace dummy {
  18   using namespace boost;
  19   using namespace numeric;
  20   using namespace interval_lib;
  21   typedef save_state<rounded_arith_opp<double> > R;
  22   typedef checking_no_nan<double, checking_no_empty<double> > P;
  23   typedef interval<double, policies<R, P> > I;
  24 }
  25
  26 template<class T>
  27 class vector {
  28   T* ptr;
  29 public:
  30   vector(int d) { ptr = (T*)malloc(sizeof(T) * d); }
  31   ~vector() { free(ptr); }
  32   const T& operator[](int i) const { return ptr[i]; }
  33   T& operator[](int i) { return ptr[i]; }
  34 };
  35
  36 template<class T>
  37 class matrix {
  38   int dim;
  39   T* ptr;
  40 public:
  41   matrix(int d): dim(d) { ptr = (T*)malloc(sizeof(T) * dim * dim); }
  42   ~matrix() { free(ptr); }
  43   int get_dim() const { return dim; }
  44   void assign(const matrix<T> &a) { memcpy(ptr, a.ptr, sizeof(T) * dim * dim); }
  45   const T* operator[](int i) const { return &(ptr[i * dim]); }
  46   T* operator[](int i) { return &(ptr[i * dim]); }
  47 };
  48
  49 typedef dummy::I I_dbl;
  50
  51 /* compute the sign of a determinant using an interval LU-decomposition; the
  52    function answers 1 or -1 if the determinant is positive or negative (and
  53    more importantly, the result must be provable), or 0 if the algorithm was
  54    unable to get a correct sign */
  55 int det_sign_algo1(const matrix<double> &a) {
  56   int dim = a.get_dim();
  57   vector<int> p(dim);
  58   for(int i = 0; i < dim; i++) p[i] = i;
  59   int sig = 1;
  60   I_dbl::traits_type::rounding rnd;
  61   typedef boost::numeric::interval_lib::unprotect<I_dbl>::type I;
  62   matrix<I> u(dim);
  63   for(int i = 0; i < dim; i++) {
  64     const double* line1 = a[i];
  65     I* line2 = u[i];
  66     for(int j = 0; j < dim; j++)
  67       line2[j] = line1[j];
  68   }
  69   // computation of L and U
  70   for(int i = 0; i < dim; i++) {
  71     // partial pivoting
  72     {
  73       int pivot = i;
  74       double max = 0;
  75       for(int j = i; j < dim; j++) {
  76         const I &v = u[p[j]][i];
  77         if (zero_in(v)) continue;
  78         double m = norm(v);
  79         if (m > max) { max = m; pivot = j; }
  80       }
  81       if (max == 0) return 0;
  82       if (pivot != i) {
  83         sig = -sig;
  84         int tmp = p[i];
  85         p[i] = p[pivot];
  86         p[pivot] = tmp;
  87       }
  88     }
  89     // U[i,?]
  90     {
  91       I *line1 = u[p[i]];
  92       const I &pivot = line1[i];
  93       if (boost::numeric::interval_lib::cerlt(pivot, 0.)) sig = -sig;
  94       for(int k = i + 1; k < dim; k++) {
  95         I *line2 = u[p[k]];
  96         I fact = line2[i] / pivot;
  97         for(int j = i + 1; j < dim; j++) line2[j] -= fact * line1[j];
  98       }
  99     }
 100   }
 101   return sig;
 102 }
 103
 104 /* compute the sign of a determinant using a floating-point LU-decomposition
 105    and an a posteriori interval validation; the meaning of the answer is the
 106    same as previously */
 107 int det_sign_algo2(const matrix<double> &a) {
 108   int dim = a.get_dim();
 109   vector<int> p(dim);
 110   for(int i = 0; i < dim; i++) p[i] = i;
 111   int sig = 1;
 112   matrix<double> lui(dim);
 113   {
 114     // computation of L and U
 115     matrix<double> lu(dim);
 116     lu.assign(a);
 117     for(int i = 0; i < dim; i++) {
 118       // partial pivoting
 119       {
 120         int pivot = i;
 121         double max = std::abs(lu[p[i]][i]);
 122         for(int j = i + 1; j < dim; j++) {
 123           double m = std::abs(lu[p[j]][i]);
 124           if (m > max) { max = m; pivot = j; }
 125         }
 126         if (max == 0) return 0;
 127         if (pivot != i) {
 128           sig = -sig;
 129           int tmp = p[i];
 130           p[i] = p[pivot];
 131           p[pivot] = tmp;
 132         }
 133       }
 134       // L[?,i] and U[i,?]
 135       {
 136         double *line1 = lu[p[i]];
 137         double pivot = line1[i];
 138         if (pivot < 0) sig = -sig;
 139         for(int k = i + 1; k < dim; k++) {
 140           double *line2 = lu[p[k]];
 141           double fact = line2[i] / pivot;
 142           line2[i] = fact;
 143           for(int j = i + 1; j < dim; j++) line2[j] -= line1[j] * fact;
 144         }
 145       }
 146     }
 147
 148     // computation of approximate inverses: Li and Ui
 149     for(int j = 0; j < dim; j++) {
 150       for(int i = j + 1; i < dim; i++) {
 151         double *line = lu[p[i]];
 152         double s = - line[j];
 153         for(int k = j + 1; k < i; k++) s -= line[k] * lui[k][j];
 154         lui[i][j] = s;
 155       }
 156       lui[j][j] = 1 / lu[p[j]][j];
 157       for(int i = j - 1; i >= 0; i--) {
 158         double *line = lu[p[i]];
 159         double s = 0;
 160         for(int k = i + 1; k <= j; k++) s -= line[k] * lui[k][j];
 161         lui[i][j] = s / line[i];
 162       }
 163     }
 164   }
 165
 166   // norm of PAUiLi-I computed with intervals
 167   {
 168     I_dbl::traits_type::rounding rnd;
 169     typedef boost::numeric::interval_lib::unprotect<I_dbl>::type I;
 170     vector<I> m1(dim);
 171     vector<I> m2(dim);
 172     for(int i = 0; i < dim; i++) {
 173       for(int j = 0; j < dim; j++) m1[j] = 0;
 174       const double *l1 = a[p[i]];
 175       for(int j = 0; j < dim; j++) {
 176         double v = l1[j];    // PA[i,j]
 177         double *l2 = lui[j]; // Ui[j,?]
 178         for(int k = j; k < dim; k++) {
 179           using boost::numeric::interval_lib::mul;
 180           m1[k] += mul<I>(v, l2[k]); // PAUi[i,k]
 181         }
 182       }
 183       for(int j = 0; j < dim; j++) m2[j] = m1[j]; // PAUi[i,j] * Li[j,j]
 184       for(int j = 1; j < dim; j++) {
 185         const I &v = m1[j];  // PAUi[i,j]
 186         double *l2 = lui[j]; // Li[j,?]
 187         for(int k = 0; k < j; k++)
 188           m2[k] += v * l2[k]; // PAUiLi[i,k]
 189       }
 190       m2[i] -= 1; // PAUiLi-I
 191       double ss = 0;
 192       for(int i = 0; i < dim; i++) ss = rnd.add_up(ss, norm(m2[i]));
 193       if (ss >= 1) return 0;
 194     }
 195   }
 196   return sig;
 197 }
 198
 199 int main() {
 200   int dim = 20;
 201   matrix<double> m(dim);
 202   for(int i = 0; i < dim; i++) for(int j = 0; j < dim; j++)
 203     m[i][j] = /*1 / (i-j-0.001)*/ cos(1+i*sin(1 + j)) /*1./(1+i+j)*/;
 204
 205   // compute the sign of the determinant of a "strange" matrix with the two
 206   // algorithms, the first should fail and the second succeed
 207   std::cout << det_sign_algo1(m) << " " << det_sign_algo2(m) << std::endl;
 208 }