Skip to content
Snippets Groups Projects
Commit d53dec64 authored by Vincent Nivoliers's avatar Vincent Nivoliers
Browse files

wrong commited file

parent cb4250de
No related branches found
No related tags found
No related merge requests found
File deleted
#include <boost/numeric/interval.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <eigen3/Eigen/Dense>
#include <random>
#include <algorithm>
#include <iostream>
//types
using interval = boost::numeric::interval<double> ;
using rational = boost::multiprecision::cpp_rational ;
using MatDX = Eigen::MatrixXd ;
using VecDX = Eigen::VectorXd ;
using MatIX = Eigen::Matrix<interval, Eigen::Dynamic, Eigen::Dynamic> ;
using VecIX = Eigen::Matrix<interval, Eigen::Dynamic, 1> ;
using MatRX = Eigen::Matrix<rational, Eigen::Dynamic, Eigen::Dynamic> ;
using VecRX = Eigen::Matrix<rational, Eigen::Dynamic, 1> ;
namespace Eigen {
namespace internal {
template<typename X, typename S, typename P>
struct is_convertible<X,boost::numeric::interval<S,P> > {
enum { value = is_convertible<X,S>::value };
};
template<typename S, typename P1, typename P2>
struct is_convertible<boost::numeric::interval<S,P1>,boost::numeric::interval<S,P2> > {
enum { value = true };
};
}
}
//command line
char* get_cmd_option( char ** begin, char ** end,
const std::string & option
) {
char ** itr = std::find(begin, end, "--" + option);
if (itr != end && ++itr != end)
{
return *itr;
}
return 0;
}
bool cmd_option_exists( char** begin, char** end,
const std::string& option
) {
return std::find(begin, end, "--" + option) != end;
}
template<typename T>
bool get_cmd_value(
char ** begin,
char** end,
const std::string& option,
T* output
) {
if(cmd_option_exists(begin, end, option)) {
std::stringstream ss ;
ss << get_cmd_option(begin, end, option) ;
T t ;
ss >> t ;
if(!ss.fail()) {
*output = t ;
return true ;
}
}
return false ;
}
int main(int argc, char** argv) {
//test variables
unsigned int dim = 10 ;
get_cmd_value(argv, argv + argc, "dim", &dim) ;
//alea
std::random_device seed ;
std::mt19937_64 mt(seed()) ;
std::uniform_real_distribution<double> random_double(-1, 1) ;
//initialization
MatDX matd(dim, dim) ;
VecDX vecd(dim) ;
MatIX mati(dim, dim) ;
VecIX veci(dim) ;
MatRX matr(dim, dim) ;
VecRX vecr(dim) ;
for(unsigned int i = 0; i < dim; ++i) {
for(unsigned int j = 0; j < dim; ++j) {
double d = random_double(mt) ;
matd(i, j) = d ;
mati(i, j) = d ;
matr(i, j) = d ;
}
double d = random_double(mt) ;
vecd(i) = d ;
veci(i) = d ;
vecr(i) = d ;
}
rational detr = matr.fullPivLu().determinant() ;
std::cout << "==========" << std::endl ;
{
std::cout << "solve fullPivLu" << std::endl ;
double detd = matd.fullPivLu().determinant() ;
try {
interval deti = mati.fullPivLu().determinant() ;
std::cout << detd
<< " in [" << deti.lower() << " -- " << deti.upper()
<< "] error " << detd - (double) detr
<< " wrt. " << deti.upper() - deti.lower() << std::endl ;
} catch(std::exception& e) {
std::cout << "interval computation failed" << std::endl ;
std::cout << detd << " error " << detd - (double) detr << std::endl ;
}
std::cout << "==========" << std::endl ;
}
{
std::cout << "solve partialPivLu" << std::endl ;
double detd = matd.partialPivLu().determinant() ;
try {
interval deti = mati.partialPivLu().determinant() ;
std::cout << detd
<< " in [" << deti.lower() << " -- " << deti.upper()
<< "] error " << detd - (double) detr
<< " wrt. " << deti.upper() - deti.lower() << std::endl ;
} catch(std::exception& e) {
std::cout << "interval computation failed" << std::endl ;
std::cout << detd << " error " << detd - (double) detr << std::endl ;
}
std::cout << "==========" << std::endl ;
}
return 0 ;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment