Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/tests/specmicp_common/condition_number.cpp
f9c5caeee1ffmaster

Last Change
Raw File

condition_number.cpp
View Options

#include "catch.hpp"
#include "specmicp_common/micpsolver/estimate_cond_number.hpp"
#include <Eigen/Dense>
#include <Eigen/Jacobi>
#include <iostream>
TEST_CASE("Condition number") {
SECTION("Lower matrix") {
Eigen::MatrixXd mat(4,4);
mat << 1, 0, 0 , 0,
1, 1, 0, 0,
1, 1, 1, 0,
1, 1, 1, 1;
double cond = specmicp::micpsolver::estimate_condition_number<Eigen::Lower>(mat);
CHECK(cond == 4);
}
SECTION("Upper matrix") {
Eigen::MatrixXd mat(4,4);
mat << 1, 1, 1 , 1,
0, 1, 1, 1,
0, 0, 1, 1,
0, 0, 0, 1;
double cond = specmicp::micpsolver::estimate_condition_number<Eigen::Upper>(mat);
CHECK(cond == 4);
}
SECTION("Lower matrix, 1000") {
Eigen::MatrixXd mat(4, 4);
mat << 1, 0, 0 , 0,
1, 1, 0, 0,
1, 1, 1, 0,
1000, 1, 1, 1;
double cond = specmicp::micpsolver::estimate_condition_number<Eigen::Lower>(mat);
CHECK(cond == 1003);
}
SECTION("Lower matrix 1000 2") {
Eigen::MatrixXd mat(4, 4);
mat << 1, 0, 0 , 0,
1, 2000, 0, 0,
1, 1, 1, 0,
1000, 1, 1, 1000;
double cond = specmicp::micpsolver::estimate_condition_number<Eigen::Lower>(mat);
CHECK(cond == 2002);
}
SECTION("Lower matrix 10000 ") {
Eigen::MatrixXd mat(4, 4);
mat << 1, 0, 0 , 0,
1, 1000, 0, 0,
1, 1, 1, 0,
1000, 1, 1, 10000;
double cond = specmicp::micpsolver::estimate_condition_number<Eigen::Lower>(mat);
CHECK(cond == 11002);
}
SECTION("Lower matrix signs ") {
Eigen::MatrixXd mat(4, 4);
mat << -1, 0, 0 , 0,
1, -2000, 0, 0,
-8, 1, -5, 0,
1000, 1, 1, -10000;
double cond = specmicp::micpsolver::estimate_condition_number<Eigen::Lower>(mat);
CHECK(cond == 11002);
}
SECTION("Random") {
int nb_mat = 30;
double sum;
for (int i=0; i<nb_mat; ++i)
{
Eigen::MatrixXd mat(30, 30);
mat.setRandom();
mat *=10;
Eigen::JacobiSVD<Eigen::MatrixXd> svd(mat);
double cond = svd.singularValues()(0)
/ svd.singularValues()(svd.singularValues().size()-1);
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr_solver(mat);
double cond_estimate = specmicp::micpsolver::estimate_condition_number<Eigen::Lower>(qr_solver.matrixR());
sum += cond_estimate/cond;
}
const double average = sum/nb_mat;
CHECK(average > 0.2);
CHECK(average < 5); // right order of magnitudes
}
}
c4science · Help