Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/tests/micpsolver/condition_number.cpp
e96f9ec8af07reactmicp/unsaturated

Last Change
Raw File

condition_number.cpp
View Options

#include "catch.hpp"
#include <Eigen/Dense>
#include "micpsolver/estimate_cond_number.hpp"
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(
mat.triangularView<Eigen::Lower>());
CHECK(cond == 4);
}
SECTION("Upper 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(
mat.transpose().triangularView<Eigen::Upper>());
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(
mat.triangularView<Eigen::Lower>());
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(
mat.triangularView<Eigen::Lower>());
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(
mat.triangularView<Eigen::Lower>());
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(
mat.triangularView<Eigen::Lower>());
CHECK(cond == 11002);
}
SECTION("Random") {
for (int i=0; i<20; ++i)
{
Eigen::MatrixXd mat(20, 20);
mat.setRandom();
mat.triangularView<Eigen::Upper>().setZero();
double norm = mat.norm();
double norm_inv = mat.inverse().norm();
double cond = specmicp::micpsolver::estimate_condition_number(
mat.triangularView<Eigen::Lower>());
CHECK(abs(cond - (norm_inv/norm)) < 1.0);
}
}
}
c4science · Help