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