Select Git revision
PugsParser.cpp
-
Stéphane Del Pino authored
The idea is to define functions this way `` let f : X -> Y, x -> y; `` where y is an expression of x and previously defined variables. - x is local to the function! - y cannot change *any* of its parameters (neither x or previously defined data) These will not be C/C++-like functions (which will be referred as routines) In this commit X and Y are basic types such as N, Z, R, B and string are allowed. The aim is to allow constructions like this, for instance: `` // definition let norm : R*R -> R, (x,y) -> sqrt(x*x + y*y); // usage R n = norm(2,3.2); `` Standard functions such as sin, cos, sqrt, abs, ... should be predefined functions of that kind. This commit only defines the simplest grammar. Definition is checked but use is still invalid (since function definition is not stored correctly according with regard to the function variable)
Stéphane Del Pino authoredThe idea is to define functions this way `` let f : X -> Y, x -> y; `` where y is an expression of x and previously defined variables. - x is local to the function! - y cannot change *any* of its parameters (neither x or previously defined data) These will not be C/C++-like functions (which will be referred as routines) In this commit X and Y are basic types such as N, Z, R, B and string are allowed. The aim is to allow constructions like this, for instance: `` // definition let norm : R*R -> R, (x,y) -> sqrt(x*x + y*y); // usage R n = norm(2,3.2); `` Standard functions such as sin, cos, sqrt, abs, ... should be predefined functions of that kind. This commit only defines the simplest grammar. Definition is checked but use is still invalid (since function definition is not stored correctly according with regard to the function variable)
test_EigenvalueSolver.cpp 5.53 KiB
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
#include <algebra/CRSMatrix.hpp>
#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/EigenvalueSolver.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
#include <scheme/HyperelasticSolver.hpp>
#include <utils/pugs_config.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("EigenvalueSolver", "[algebra]")
{
SECTION("Sparse Matrices")
{
SECTION("symmetric system")
{
Array<int> non_zeros(3);
non_zeros.fill(3);
CRSMatrixDescriptor<double> S{3, 3, non_zeros};
S(0, 0) = 3;
S(0, 1) = 2;
S(0, 2) = 4;
S(1, 0) = 2;
S(1, 1) = 0;
S(1, 2) = 2;
S(2, 0) = 4;
S(2, 1) = 2;
S(2, 2) = 3;
CRSMatrix A{S.getCRSMatrix()};
SECTION("eigenvalues")
{
SmallArray<double> eigenvalues;
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues);
REQUIRE(eigenvalues[0] == Catch::Approx(-1));
REQUIRE(eigenvalues[1] == Catch::Approx(-1));
REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else // PUGS_HAS_SLEPC
REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues),
"not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
SECTION("eigenvalues and eigenvectors")
{
SmallArray<double> eigenvalue_list;
std::vector<SmallVector<double>> eigenvector_list;
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
SmallMatrix<double> P{3};
SmallMatrix<double> PT{3};
SmallMatrix<double> D{3};
D = zero;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
P(i, j) = eigenvector_list[j][i];
PT(i, j) = eigenvector_list[i][j];
}
D(i, i) = eigenvalue_list[i];
}
SmallMatrix<double> eigen{3};
eigen = zero;
SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
}
}
#else // PUGS_HAS_SLEPC
REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
"not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
SECTION("eigenvalues and passage matrix")
{
SmallArray<double> eigenvalue_list;
SmallMatrix<double> P{};
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);
REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
SmallMatrix<double> D{3};
D = zero;
for (size_t i = 0; i < 3; ++i) {
D(i, i) = eigenvalue_list[i];
}
SmallMatrix PT = transpose(P);
TinyMatrix<3, 3, double> TinyA;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
TinyA(i, j) = S(i, j);
}
}
SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
}
}
#else // PUGS_HAS_SLEPC
REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P),
"not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
SECTION("exponential of a matrix")
{
SmallMatrix<double> expA{};
SmallMatrix<double> expS{3};
expS(0, 0) = 1325.074593930307812;
expS(0, 1) = 662.353357244568285;
expS(0, 2) = 1324.70671448913637;
expS(1, 0) = expS(0, 1);
expS(1, 1) = 331.544558063455535;
expS(1, 2) = 662.353357244568185; expS(2, 0) = expS(0, 2);
expS(2, 1) = expS(1, 2);
expS(2, 2) = expS(0, 0);
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA);
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
REQUIRE(expA(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12));
}
}
#else // PUGS_HAS_SLEPC
REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA),
"not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif // PUGS_HAS_SLEPC
}
}
SECTION("symmetric tiny matrix")
{
TinyMatrix<3> TestA{3e10, 2e10, 4e10, 2e10, 0, 2e10, 4e10, 2e10, 3e10};
TinyMatrix<3> TestA2{4, 2, 3, 2, 0, 2, 3, 2, 4};
TinyMatrix<3> TestB{1, -1, 0, -1, 1, 0, 0, 0, 3};
TinyMatrix<3> TestC{0, 0, 0, 0, 0, 0, 0, 0, 0};
TinyMatrix<3> TestD{1e10, 0, 0, 0, 1, 0, 0, 0, 1};
TinyMatrix<3> TestE{3e-10, 2e-10, 4e-10, 2e-10, 0, 2e-10, 4e-10, 2e-10, 3e-10};
TinyMatrix<3> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1};
TinyVector<3> eigenvalues0{0, 0, 0};
auto [eigenvalues, eigenmatrix] = EigenvalueSolver{}.findEigen(TestA);
}
}
}