interval/examples/filter.cpp

/* 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 <cstring>
#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;
}