GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-11 03:09:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

added functions to fetch L and U Factors from IncompleteLUT

Browse Source
This commit is contained in:
Johannes Zipfel 2025-01-31 18:32:38 +00:00 committed by Rasmus Munk Larsen
parent b6849f675d
commit 2926b2e0a9
3 changed files with 119 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
View File

@ -129,6 +129,12 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<Scalar_, StorageInde
compute(mat);
}
/** \brief Extraction Method for L-Factor */
const FactorType matrixL() const;
/** \brief Extraction Method for U-Factor */
const FactorType matrixU() const;
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
@ -207,6 +213,28 @@ void IncompleteLUT<Scalar, StorageIndex>::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 <typename Scalar, typename StorageIndex>
const typename IncompleteLUT<Scalar, StorageIndex>::FactorType IncompleteLUT<Scalar, StorageIndex>::matrixL() const {
eigen_assert(m_factorizationIsOk && "factorize() should be called first");
return m_lu.template triangularView<UnitLower>();
}
/**
* 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 <typename Scalar, typename StorageIndex>
const typename IncompleteLUT<Scalar, StorageIndex>::FactorType IncompleteLUT<Scalar, StorageIndex>::matrixU() const {
eigen_assert(m_factorizationIsOk && "Factorization must be computed first.");
return m_lu.template triangularView<Upper>();
}
template <typename Scalar, typename StorageIndex>
template <typename MatrixType_>
void IncompleteLUT<Scalar, StorageIndex>::analyzePattern(const MatrixType_& amat) {
@ -418,4 +446,4 @@ void IncompleteLUT<Scalar, StorageIndex>::factorize(const MatrixType_& amat) {
} // end namespace Eigen
#endif // EIGEN_INCOMPLETE_LUT_H
#endif // EIGEN_INCOMPLETE_LUT_H

1
test/CMakeLists.txt
View File

@ -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)

89
test/incomplete_LUT.cpp Normal file
View File

@ -0,0 +1,89 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2025 Johannes Zipfel <johzip1010@gmail.com>
//
// 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 <Eigen/IterativeLinearSolvers>
template <typename T, typename I_>
void test_incompleteLUT_T() {
IncompleteLUT<T, I_> ilut;
ilut.setDroptol(NumTraits<T>::epsilon() * 4);
}
template <typename T>
void test_extract_LU() {
typedef Eigen::SparseMatrix<T> SparseMatrix;
SparseMatrix A(5, 5);
std::vector<Eigen::Triplet<T>> 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<T> ilut;
ilut.compute(A);
Eigen::SparseMatrix<T> matL = ilut.matrixL(); // Extract L
Eigen::SparseMatrix<T> matU = ilut.matrixU(); // Extract U
Eigen::SparseMatrix<T> expectedMatL(5, 5);
std::vector<Eigen::Triplet<T>> 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<T> expectedMatU(5, 5);
std::vector<Eigen::Triplet<T>> 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<double, int>()));
CALL_SUBTEST_1((test_incompleteLUT_T<float, int>()));
CALL_SUBTEST_2((test_incompleteLUT_T<std::complex<double>, int>()));
CALL_SUBTEST_3((test_incompleteLUT_T<double, long int>()));
CALL_SUBTEST_4(test_extract_LU<double>());
CALL_SUBTEST_4(test_extract_LU<float>());
}