From 2926b2e0a91dcd4ca9b01758a12b73ccfc91894b Mon Sep 17 00:00:00 2001 From: Johannes Zipfel Date: Fri, 31 Jan 2025 18:32:38 +0000 Subject: [PATCH] added functions to fetch L and U Factors from IncompleteLUT --- .../IterativeLinearSolvers/IncompleteLUT.h | 30 ++++++- test/CMakeLists.txt | 1 + test/incomplete_LUT.cpp | 89 +++++++++++++++++++ 3 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 test/incomplete_LUT.cpp diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 575a7b223..930077df4 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -129,6 +129,12 @@ class IncompleteLUT : public SparseSolverBase::setFillfactor(int fillfactor) { this->m_fillfactor = fillfactor; } +/** + * get L-Factor + * \return L-Factor is a matrix containing the lower triangular part of the sparse matrix. All elements of the matrix + * above the main diagonal are zero. + **/ +template +const typename IncompleteLUT::FactorType IncompleteLUT::matrixL() const { + eigen_assert(m_factorizationIsOk && "factorize() should be called first"); + return m_lu.template triangularView(); +} + +/** + * get U-Factor + * \return L-Factor is a matrix containing the upper triangular part of the sparse matrix. All elements of the matrix + * below the main diagonal are zero. + **/ +template +const typename IncompleteLUT::FactorType IncompleteLUT::matrixU() const { + eigen_assert(m_factorizationIsOk && "Factorization must be computed first."); + return m_lu.template triangularView(); +} + template template void IncompleteLUT::analyzePattern(const MatrixType_& amat) { @@ -418,4 +446,4 @@ void IncompleteLUT::factorize(const MatrixType_& amat) { } // end namespace Eigen -#endif // EIGEN_INCOMPLETE_LUT_H +#endif // EIGEN_INCOMPLETE_LUT_H \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 469258435..fdfde4597 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -285,6 +285,7 @@ ei_add_test(sparse_permutations) ei_add_test(simplicial_cholesky) ei_add_test(conjugate_gradient) ei_add_test(incomplete_cholesky) +ei_add_test(incomplete_LUT) ei_add_test(bicgstab) ei_add_test(lscg) ei_add_test(sparselu) diff --git a/test/incomplete_LUT.cpp b/test/incomplete_LUT.cpp new file mode 100644 index 000000000..7213eccce --- /dev/null +++ b/test/incomplete_LUT.cpp @@ -0,0 +1,89 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2025 Johannes Zipfel +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "sparse_solver.h" +#include + +template +void test_incompleteLUT_T() { + IncompleteLUT ilut; + ilut.setDroptol(NumTraits::epsilon() * 4); +} + +template +void test_extract_LU() { + typedef Eigen::SparseMatrix SparseMatrix; + + SparseMatrix A(5, 5); + std::vector> triplets; + triplets.push_back({0, 0, 4}); + triplets.push_back({0, 1, -1}); + triplets.push_back({0, 4, -1}); + triplets.push_back({1, 0, -1}); + triplets.push_back({1, 1, 4}); + triplets.push_back({1, 2, -1}); + triplets.push_back({2, 1, -1}); + triplets.push_back({2, 2, 4}); + triplets.push_back({2, 3, -1}); + triplets.push_back({3, 2, -1}); + triplets.push_back({3, 3, 4}); + triplets.push_back({3, 4, -1}); + triplets.push_back({4, 0, -1}); + triplets.push_back({4, 3, -1}); + triplets.push_back({4, 4, 4}); + + A.setFromTriplets(triplets.begin(), triplets.end()); + + IncompleteLUT ilut; + ilut.compute(A); + + Eigen::SparseMatrix matL = ilut.matrixL(); // Extract L + Eigen::SparseMatrix matU = ilut.matrixU(); // Extract U + + Eigen::SparseMatrix expectedMatL(5, 5); + std::vector> tripletsExL; + tripletsExL.emplace_back(0, 0, 1); + tripletsExL.emplace_back(1, 0, -0.25); + tripletsExL.emplace_back(1, 1, 1); + tripletsExL.emplace_back(2, 0, -0.25); + tripletsExL.emplace_back(2, 1, -0.0666667); + tripletsExL.emplace_back(2, 2, 1); + tripletsExL.emplace_back(3, 2, -0.25); + tripletsExL.emplace_back(3, 3, 1); + tripletsExL.emplace_back(4, 1, -0.266667); + tripletsExL.emplace_back(4, 3, -0.266667); + tripletsExL.emplace_back(4, 4, 1); + expectedMatL.setFromTriplets(tripletsExL.begin(), tripletsExL.end()); + + Eigen::SparseMatrix expectedMatU(5, 5); + std::vector> tripletsExU; + tripletsExU.emplace_back(0, 0, 4); + tripletsExU.emplace_back(0, 1, -1); + tripletsExU.emplace_back(1, 1, 3.75); + tripletsExU.emplace_back(1, 4, -1); + tripletsExU.emplace_back(2, 2, 4); + tripletsExU.emplace_back(2, 3, -1); + tripletsExU.emplace_back(3, 3, 3.75); + tripletsExU.emplace_back(3, 4, -1); + tripletsExU.emplace_back(4, 4, 3.46667); + expectedMatU.setFromTriplets(tripletsExU.begin(), tripletsExU.end()); + + VERIFY_IS_APPROX(expectedMatL, matL); + VERIFY_IS_APPROX(expectedMatU, matU); +} + +EIGEN_DECLARE_TEST(incomplete_LUT) { + CALL_SUBTEST_1((test_incompleteLUT_T())); + CALL_SUBTEST_1((test_incompleteLUT_T())); + CALL_SUBTEST_2((test_incompleteLUT_T, int>())); + CALL_SUBTEST_3((test_incompleteLUT_T())); + + CALL_SUBTEST_4(test_extract_LU()); + CALL_SUBTEST_4(test_extract_LU()); +}