Skip to content

L
libcifpp
Commit 76c5706f by Maarten L. Hekkelman
Options

moving to eigen3 eigensolver, fixing include and dependencies

parent 2bf4284f
Showing with 4 additions and 184 deletions
CMakeLists.txt
...@@ -168,7 +168,7 @@ if(MSVC) ...@@ -168,7 +168,7 @@ if(MSVC)
endif() endif()
find_package(ZLIB REQUIRED) find_package(ZLIB REQUIRED)
find_package(Eigen3 ) find_package(Eigen3 REQUIRED)
include(FindFilesystem) include(FindFilesystem)
list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY}) list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY})
...@@ -275,11 +275,11 @@ set_target_properties(cifpp PROPERTIES POSITION_INDEPENDENT_CODE ON) ...@@ -275,11 +275,11 @@ set_target_properties(cifpp PROPERTIES POSITION_INDEPENDENT_CODE ON)
target_include_directories(cifpp target_include_directories(cifpp
PUBLIC PUBLIC
"$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include;${PROJECT_BINARY_DIR}>" "$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include;${PROJECT_BINARY_DIR};${EIGEN3_INCLUDE_DIR}>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>" "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
) )
target_link_libraries(cifpp PUBLIC Threads::Threads ZLIB::ZLIB ${CIFPP_REQUIRED_LIBRARIES} Eigen3::Eigen) target_link_libraries(cifpp PUBLIC Threads::Threads ZLIB::ZLIB ${CIFPP_REQUIRED_LIBRARIES})
if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang") if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
target_link_options(cifpp PRIVATE -undefined dynamic_lookup) target_link_options(cifpp PRIVATE -undefined dynamic_lookup)
...@@ -426,7 +426,7 @@ if(ENABLE_TESTING) ...@@ -426,7 +426,7 @@ if(ENABLE_TESTING)
add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE}) add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE})
target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp::cifpp Boost::boost Eigen3::Eigen) target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp::cifpp Boost::boost)
if(MSVC) if(MSVC)
# Specify unwind semantics so that MSVC knowns how to handle exceptions # Specify unwind semantics so that MSVC knowns how to handle exceptions
......
include/cif++/matrix.hpp
...@@ -473,137 +473,6 @@ matrix3x3<F> inverse(const matrix3x3<F> &m) ...@@ -473,137 +473,6 @@ matrix3x3<F> inverse(const matrix3x3<F> &m)
return result; return result;
} }
template <typename M>
auto eigen(const matrix_expression<M> &mat, bool sort)
{
using value_type = decltype(mat(0, 0));
assert(mat.dim_m() == mat.dim_n());
const size_t N = mat.dim_m();
matrix<float> m = mat;
matrix<value_type> em = identity_matrix(N);
std::vector<value_type> ev(N);
std::vector<value_type> b(N);
// Set ev & b to diagonal.
for (size_t i = 0; i < N; ++i)
ev[i] = b[i] = m(i, i);
for (int cyc = 1; cyc <= 50; ++cyc)
{
// calc sum of diagonal, off-diagonal
value_type spp = 0, spq = 0;
for (size_t i = 0; i < N - 1; i++)
{
for (size_t j = i + 1; j < N; j++)
spq += std::fabs(m(i, j));
spp += std::fabs(m(i, i));
}
if (spq <= 1.0e-12 * spp)
break;
std::vector<value_type> z(N);
// now try and reduce each off-diagonal element in turn
for (size_t i = 0; i < N - 1; i++)
{
for (size_t j = i + 1; j < N; j++)
{
value_type a_ij = m(i, j);
value_type h = ev[j] - ev[i];
value_type t;
if (h == 0 and a_ij == 0)
t = 1;
else if (std::fabs(a_ij) > 1.0e-12 * std::fabs(h))
{
value_type theta = 0.5 * h / a_ij;
t = 1.0 / (std::fabs(theta) + std::sqrt(1 + theta * theta));
if (theta < 0)
t = -t;
}
else
t = a_ij / h;
// calc trig properties
value_type c = 1.f / std::sqrt(1 + t * t);
value_type s = t * c;
value_type tau = s / (1 + c);
h = t * a_ij;
// update eigenvalues
z[i] -= h;
z[j] += h;
ev[i] -= h;
ev[j] += h;
// rotate the upper diagonal of the matrix
m(i, j) = 0;
for (size_t k = 0; k < i; k++)
{
value_type ai = m(k, i);
value_type aj = m(k, j);
m(k, i) = ai - s * (aj + ai * tau);
m(k, j) = aj + s * (ai - aj * tau);
}
for (size_t k = i + 1; k < j; k++)
{
value_type ai = m(i, k);
value_type aj = m(k, j);
m(i, k) = ai - s * (aj + ai * tau);
m(k, j) = aj + s * (ai - aj * tau);
}
for (size_t k = j + 1; k < N; k++)
{
value_type ai = m(i, k);
value_type aj = m(j, k);
m(i, k) = ai - s * (aj + ai * tau);
m(j, k) = aj + s * (ai - aj * tau);
}
// apply corresponding rotation to result
for (size_t k = 0; k < N; k++)
{
value_type ai = em(k, i);
value_type aj = em(k, j);
em(k, i) = ai - s * (aj + ai * tau);
em(k, j) = aj + s * (ai - aj * tau);
}
}
}
for (size_t p = 0; p < N; p++)
{
b[p] += z[p];
ev[p] = b[p];
}
}
if (sort)
{
for (size_t p = 0; p < N; ++p)
{
size_t j = p; // set j to index of largest remaining eval
for (size_t q = p + 1; q < N; ++q)
if (ev[q] < ev[j])
j = q;
std::swap(ev[j], ev[p]);
em.swap_col(j, p);
}
}
return std::make_tuple(ev, em);
}
// -------------------------------------------------------------------- // --------------------------------------------------------------------
template <typename M> template <typename M>
......
test/unit-3d-test.cpp
...@@ -481,52 +481,4 @@ BOOST_AUTO_TEST_CASE(symm_3bwh_1, *utf::tolerance(0.1f)) ...@@ -481,52 +481,4 @@ BOOST_AUTO_TEST_CASE(symm_3bwh_1, *utf::tolerance(0.1f))
BOOST_TEST(d == distance(a1.get_location(), p)); BOOST_TEST(d == distance(a1.get_location(), p));
} }
} }
}
// --------------------------------------------------------------------
BOOST_AUTO_TEST_CASE(eigen_1, *utf::tolerance(0.1f))
{
cif::symmetric_matrix4x4<float> m;
m(0, 0) = 4;
m(0, 1) = -30;
m(0, 2) = 60;
m(0, 3) = -35;
m(1, 1) = 300;
m(1, 2) = -675;
m(1, 3) = 420;
m(2, 2) = 1620;
m(2, 3) = -1050;
m(3, 3) = 700;
cif::matrix4x4<float> m2;
m2 = m;
const auto &[ev, em] = cif::eigen(m2, true);
BOOST_TEST(ev[0] == 0.1666428611718905f);
BOOST_TEST(ev[1] == 1.4780548447781369f);
BOOST_TEST(ev[2] == 37.1014913651276582f);
BOOST_TEST(ev[3] == 2585.25381092892231f);
BOOST_TEST(em(0, 0) == 0.792608291163763585f);
BOOST_TEST(em(1, 0) == 0.451923120901599794f);
BOOST_TEST(em(2, 0) == 0.322416398581824992f);
BOOST_TEST(em(3, 0) == 0.252161169688241933f);
BOOST_TEST(em(0, 1) == -0.582075699497237650f);
BOOST_TEST(em(1, 1) == 0.370502185067093058f);
BOOST_TEST(em(2, 1) == 0.509578634501799626f);
BOOST_TEST(em(3, 1) == 0.514048272222164294f);
// BOOST_TEST(em(0, 2) == -0.179186290535454826f);
// BOOST_TEST(em(1, 2) == 0.741917790628453435f);
// BOOST_TEST(em(2, 2) == -0.100228136947192199f);
// BOOST_TEST(em(3, 2) == -0.638282528193614892f);
BOOST_TEST(em(0, 3) == 0.0291933231647860588f);
BOOST_TEST(em(1, 3) == -0.328712055763188997f);
BOOST_TEST(em(2, 3) == 0.791411145833126331f);
BOOST_TEST(em(3, 3) == -0.514552749997152907f);
} }
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment