From 618df91c5af399c26772e0418fdc95b6942e125d Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 9 Jan 2025 15:14:33 -0800 Subject: [PATCH 1/2] Add in-place LU decomposition with Doolittle algorithm (#705) * add in-place Doolittle LU decomposition * add CPU solver builder for in-place LU decomp * add typename --- include/micm/solver/lu_decomposition.hpp | 9 + .../lu_decomposition_doolittle_in_place.hpp | 106 +++++++++ .../lu_decomposition_doolittle_in_place.inl | 220 ++++++++++++++++++ include/micm/solver/solver_builder.hpp | 23 ++ .../test_analytical_backward_euler.cpp | 70 ++---- test/unit/solver/CMakeLists.txt | 1 + ...st_lu_decomposition_doolittle_in_place.cpp | 75 ++++++ 7 files changed, 458 insertions(+), 46 deletions(-) create mode 100644 include/micm/solver/lu_decomposition_doolittle_in_place.hpp create mode 100644 include/micm/solver/lu_decomposition_doolittle_in_place.inl create mode 100644 test/unit/solver/test_lu_decomposition_doolittle_in_place.cpp diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 5bcdcb1bc..9c04fb8c8 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -4,6 +4,7 @@ #include "lu_decomposition_doolittle.hpp" #include "lu_decomposition_mozart.hpp" +#include "lu_decomposition_doolittle_in_place.hpp" #include "lu_decomposition_mozart_in_place.hpp" #include @@ -26,6 +27,14 @@ namespace micm LuDecompositionMozartInPlace, SparseMatrix>>, "LuDecompositionMozartInPlace for vector matrices does not meet the LuDecompositionInPlaceConcept requirements"); + static_assert( + LuDecompositionInPlaceConcept, + "LuDecompositionDoolittleInPlace does not meet the LuDecompositionInPlaceConcept requirements"); + static_assert( + LuDecompositionInPlaceConcept< + LuDecompositionDoolittleInPlace, + SparseMatrix>>, + "LuDecompositionDoolittleInPlace for vector matrices does not meet the LuDecompositionInPlaceConcept requirements"); /// @brief Alias for the default LU decomposition algorithm using LuDecomposition = LuDecompositionDoolittle; diff --git a/include/micm/solver/lu_decomposition_doolittle_in_place.hpp b/include/micm/solver/lu_decomposition_doolittle_in_place.hpp new file mode 100644 index 000000000..7c801f295 --- /dev/null +++ b/include/micm/solver/lu_decomposition_doolittle_in_place.hpp @@ -0,0 +1,106 @@ +// Copyright (C) 2023-2024 National Center for Atmospheric Research +// SPDX-License-Identifier: Apache-2.0 +#pragma once + +#include +#include + +namespace micm +{ + + /// @brief LU decomposer for SparseMatrix following the Doolittle algorithm + /// + /// The LU decomposition uses the Doolittle algorithm following the + /// naming used here: https://www.geeksforgeeks.org/doolittle-algorithm-lu-decomposition/ + /// + /// The sudo-code for the corresponding dense matrix algorithm for matrix A (in-line) would be: + /// + /// for i = 0...n-1 // Outer loop over rows (columns) for upper (lower) triangular matrix + /// for k = i...n-1 // Middle loop over columns for upper triangular matrix + /// for j = 0...i-1 // Inner loop over columns (rows) for lower (upper) triangular matrix + /// A[i][k] -= A[i][j] * A[j][k] + /// for k = i+1...n-1 // Middle loop over rows for lower triangular matrix + /// for j = 0...i-1 // Inner loop over columns (rows) for lower (upper) triangular matrix + /// A[k][i] -= A[k][j] * A[j][i]; + /// A[k][i] /= A[i][i] + /// + /// For the sparse matrix algorithm, the indices of non-zero terms are stored in + /// several arrays during construction. These arrays are iterated through during + /// calls to Decompose to do the actual decomposition. + /// Our LU Decomposition only assigns the values of the jacobian to the LU matrices + /// when the *jacobian* is nonzero. However, the sparsity pattern of the jacobian doesn't + /// necessarily match that of the LU matrices. There can be more nonzero elements in the LU matrices + /// than in the jacobian. It is expected that the elements of the L and U matrices that are zero in the A matrix + /// will be set to zero before the combined matrix is passed to the decomposition function. + class LuDecompositionDoolittleInPlace + { + protected: + /// Number of elements in the middle (k) loops for lower and upper triangular matrices, respectively, + /// and the index in A.data_ for A[i][i] for each iteration of the outer (i) loop + std::vector> nik_nki_aii_; + /// Index in A.data_ for A[i][k] for each iteration of the upper middle (k) loop and the + /// number of elements in the inner (j) loop for each upper (k) element used to set A[i][k] + std::vector> aik_njk_; + /// Index in A.data_ for A[i][j] and A[j][k] for each iteration of the upper inner (j) loop + std::vector> aij_ajk_; + /// Index in A.data_ for A[k][i] for each iteration of the lower middle (k) loop, the + /// number of elements in the inner (j) loop for each lower (k) element used to set A[k][i] + std::vector> aki_nji_; + /// Index in A.data_ for A[k][j] and A[j][i] for each iteration of the lower inner (j) loop + std::vector> akj_aji_; + + public: + /// @brief default constructor + LuDecompositionDoolittleInPlace(); + + LuDecompositionDoolittleInPlace(const LuDecompositionDoolittleInPlace&) = delete; + LuDecompositionDoolittleInPlace& operator=(const LuDecompositionDoolittleInPlace&) = delete; + + LuDecompositionDoolittleInPlace(LuDecompositionDoolittleInPlace&& other) = default; + LuDecompositionDoolittleInPlace& operator=(LuDecompositionDoolittleInPlace&&) = default; + + /// @brief Construct an LU decomposition algorithm for a given sparse matrix + /// @param matrix Sparse matrix + template + requires(SparseMatrixConcept) + LuDecompositionDoolittleInPlace(const SparseMatrixPolicy& matrix); + + ~LuDecompositionDoolittleInPlace() = default; + + /// @brief Create an LU decomposition algorithm for a given sparse matrix policy + /// @param matrix Sparse matrix + template + requires(SparseMatrixConcept) + static LuDecompositionDoolittleInPlace Create(const SparseMatrixPolicy& matrix); + + /// @brief Create sparse L and U matrices for a given A matrix + /// @param A Sparse matrix that will be decomposed + /// @return L and U Sparse matrices + template + requires(SparseMatrixConcept) + static SparseMatrixPolicy GetLUMatrix( + const SparseMatrixPolicy& A, + typename SparseMatrixPolicy::value_type initial_value); + + /// @brief Perform an LU decomposition on a given A matrix + /// @param A Sparse matrix to decompose + /// @param L The lower triangular matrix created by decomposition + /// @param U The upper triangular matrix created by decomposition + template + requires(!VectorizableSparse) + void Decompose(SparseMatrixPolicy& ALU) const; + template + requires(VectorizableSparse) + void Decompose(SparseMatrixPolicy& ALU) const; + + protected: + /// @brief Initialize arrays for the LU decomposition + /// @param A Sparse matrix to decompose + template + requires(SparseMatrixConcept) + void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); + }; + +} // namespace micm + +#include "lu_decomposition_doolittle_in_place.inl" diff --git a/include/micm/solver/lu_decomposition_doolittle_in_place.inl b/include/micm/solver/lu_decomposition_doolittle_in_place.inl new file mode 100644 index 000000000..514ead06c --- /dev/null +++ b/include/micm/solver/lu_decomposition_doolittle_in_place.inl @@ -0,0 +1,220 @@ +// Copyright (C) 2023-2024 National Center for Atmospheric Research +// SPDX-License-Identifier: Apache-2.0 + +namespace micm +{ + + inline LuDecompositionDoolittleInPlace::LuDecompositionDoolittleInPlace() + { + } + + template + requires(SparseMatrixConcept) + inline LuDecompositionDoolittleInPlace::LuDecompositionDoolittleInPlace(const SparseMatrixPolicy& matrix) + { + Initialize(matrix, typename SparseMatrixPolicy::value_type()); + } + + template + requires(SparseMatrixConcept) + inline LuDecompositionDoolittleInPlace LuDecompositionDoolittleInPlace::Create(const SparseMatrixPolicy& matrix) + { + LuDecompositionDoolittleInPlace lu_decomp{}; + lu_decomp.Initialize(matrix, typename SparseMatrixPolicy::value_type()); + return lu_decomp; + } + + template + requires(SparseMatrixConcept) + inline void LuDecompositionDoolittleInPlace::Initialize(const SparseMatrixPolicy& matrix, auto initial_value) + { + MICM_PROFILE_FUNCTION(); + + std::size_t n = matrix.NumRows(); + auto ALU = GetLUMatrix(matrix, initial_value); + for (std::size_t i = 0; i < n; ++i) + { + if (ALU.IsZero(i, i)) + { + throw std::runtime_error("Diagonal element is zero in LU decomposition"); + } + std::tuple nik_nki_aii(0, 0, ALU.VectorIndex(0, i, i)); + for (std::size_t k = i; k < n; ++k) + { + if (ALU.IsZero(i, k)) + continue; + std::pair aik_njk(ALU.VectorIndex(0, i, k), 0); + for (std::size_t j = 0; j < i; ++j) + { + if (ALU.IsZero(i, j) || ALU.IsZero(j, k)) + continue; + aij_ajk_.push_back(std::make_pair(ALU.VectorIndex(0, i, j), ALU.VectorIndex(0, j, k))); + ++(aik_njk.second); + } + aik_njk_.push_back(aik_njk); + ++(std::get<0>(nik_nki_aii)); + } + for (std::size_t k = i + 1; k < n; ++k) + { + if (ALU.IsZero(k, i)) + continue; + std::pair aki_nji(ALU.VectorIndex(0, k, i), 0); + for (std::size_t j = 0; j < i; ++j) + { + if (ALU.IsZero(k, j) || ALU.IsZero(j, i)) + continue; + akj_aji_.push_back(std::make_pair(ALU.VectorIndex(0, k, j), ALU.VectorIndex(0, j, i))); + ++(aki_nji.second); + } + aki_nji_.push_back(aki_nji); + ++(std::get<1>(nik_nki_aii)); + } + nik_nki_aii_.push_back(nik_nki_aii); + } + } + + template + requires(SparseMatrixConcept) + inline SparseMatrixPolicy LuDecompositionDoolittleInPlace::GetLUMatrix( + const SparseMatrixPolicy& A, + typename SparseMatrixPolicy::value_type initial_value) + { + MICM_PROFILE_FUNCTION(); + + std::size_t n = A.NumRows(); + std::set> ALU_ids; + for (std::size_t i = 0; i < n; ++i) + { + // Upper triangular matrix + for (std::size_t k = i; k < n; ++k) + { + if (!A.IsZero(i, k) || k == i) + { + ALU_ids.insert(std::make_pair(i, k)); + continue; + } + for (std::size_t j = 0; j < i; ++j) + { + if (ALU_ids.find(std::make_pair(i, j)) != ALU_ids.end() && ALU_ids.find(std::make_pair(j, k)) != ALU_ids.end()) + { + ALU_ids.insert(std::make_pair(i, k)); + break; + } + } + } + // Lower triangular matrix + for (std::size_t k = i; k < n; ++k) + { + if (!A.IsZero(k, i) || k == i) + { + ALU_ids.insert(std::make_pair(k, i)); + continue; + } + for (std::size_t j = 0; j < i; ++j) + { + if (ALU_ids.find(std::make_pair(k, j)) != ALU_ids.end() && ALU_ids.find(std::make_pair(j, i)) != ALU_ids.end()) + { + ALU_ids.insert(std::make_pair(k, i)); + break; + } + } + } + } + auto ALU_builder = SparseMatrixPolicy::Create(n).SetNumberOfBlocks(A.NumberOfBlocks()).InitialValue(initial_value); + for (auto& pair : ALU_ids) + { + ALU_builder = ALU_builder.WithElement(pair.first, pair.second); + } + return SparseMatrixPolicy(ALU_builder); + } + + template + requires(!VectorizableSparse) + inline void LuDecompositionDoolittleInPlace::Decompose(SparseMatrixPolicy& ALU) const + { + MICM_PROFILE_FUNCTION(); + const std::size_t n = ALU.NumRows(); + + // Loop over blocks + for (std::size_t i_block = 0; i_block < ALU.NumberOfBlocks(); ++i_block) + { + auto ALU_vector = std::next(ALU.AsVector().begin(), i_block * ALU.FlatBlockSize()); + auto aik_njk = aik_njk_.begin(); + auto aij_ajk = aij_ajk_.begin(); + auto aki_nji = aki_nji_.begin(); + auto akj_aji = akj_aji_.begin(); + + for (const auto& nik_nki_aii : nik_nki_aii_) + { + for (std::size_t ik = 0; ik < std::get<0>(nik_nki_aii); ++ik) + { + for (std::size_t jk = 0; jk < aik_njk->second; ++jk) + { + ALU_vector[aik_njk->first] -= ALU_vector[aij_ajk->first] * ALU_vector[aij_ajk->second]; + ++aij_ajk; + } + ++aik_njk; + } + for (std::size_t ki = 0; ki < std::get<1>(nik_nki_aii); ++ki) + { + for (std::size_t ji = 0; ji < aki_nji->second; ++ji) + { + ALU_vector[aki_nji->first] -= ALU_vector[akj_aji->first] * ALU_vector[akj_aji->second]; + ++akj_aji; + } + ALU_vector[aki_nji->first] /= ALU_vector[std::get<2>(nik_nki_aii)]; + ++aki_nji; + } + } + } + } + + template + requires(VectorizableSparse) + inline void LuDecompositionDoolittleInPlace::Decompose(SparseMatrixPolicy& ALU) const + { + MICM_PROFILE_FUNCTION(); + + const std::size_t n = ALU.NumRows(); + const std::size_t ALU_BlockSize = ALU.NumberOfBlocks(); + constexpr std::size_t ALU_GroupVectorSize = SparseMatrixPolicy::GroupVectorSize(); + const std::size_t ALU_GroupSizeOfFlatBlockSize = ALU.GroupSize(); + + // Loop over groups of blocks + for (std::size_t i_group = 0; i_group < ALU.NumberOfGroups(ALU_BlockSize); ++i_group) + { + auto ALU_vector = std::next(ALU.AsVector().begin(), i_group * ALU_GroupSizeOfFlatBlockSize); + const std::size_t n_cells = std::min(ALU_GroupVectorSize, ALU_BlockSize - i_group * ALU_GroupVectorSize); + auto aik_njk = aik_njk_.begin(); + auto aij_ajk = aij_ajk_.begin(); + auto aki_nji = aki_nji_.begin(); + auto akj_aji = akj_aji_.begin(); + for (const auto& nik_nki_aii : nik_nki_aii_) + { + for (std::size_t ik = 0; ik < std::get<0>(nik_nki_aii); ++ik) + { + for (std::size_t jk = 0; jk < aik_njk->second; ++jk) + { + for (std::size_t i = 0; i < n_cells; ++i) + ALU_vector[aik_njk->first + i] -= ALU_vector[aij_ajk->first + i] * ALU_vector[aij_ajk->second + i]; + ++aij_ajk; + } + ++aik_njk; + } + for (std::size_t ki = 0; ki < std::get<1>(nik_nki_aii); ++ki) + { + for (std::size_t ji = 0; ji < aki_nji->second; ++ji) + { + for (std::size_t i = 0; i < n_cells; ++i) + ALU_vector[aki_nji->first + i] -= ALU_vector[akj_aji->first + i] * ALU_vector[akj_aji->second + i]; + ++akj_aji; + } + for (std::size_t i = 0; i < n_cells; ++i) + ALU_vector[aki_nji->first + i] /= ALU_vector[std::get<2>(nik_nki_aii) + i]; + ++aki_nji; + } + } + } + } + +} // namespace micm diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index 573f4a51c..0f6c33f2f 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -42,6 +42,8 @@ namespace micm using DenseMatrixPolicyType = DenseMatrixPolicy; using SparseMatrixPolicyType = SparseMatrixPolicy; using LuDecompositionPolicyType = LuDecompositionPolicy; + using LinearSolverPolicyType = LinearSolverPolicy; + using StatePolicyType = StatePolicy; protected: SolverParametersPolicy options_; @@ -116,6 +118,8 @@ namespace micm /// @tparam DenseMatrixPolicy Policy for dense matrices /// @tparam SparseMatrixPolicy Policy for sparse matrices /// @tparam LuDecompositionPolicy Policy for the LU decomposition + /// @tparam LMatrixPolicy Policy for the Lower matrix + /// @tparam UMatrixPolicy Policy for the Upper matrix template< class SolverParametersPolicy, class DenseMatrixPolicy = Matrix, @@ -132,6 +136,25 @@ namespace micm LinearSolver, State>; + /// @brief Builder of CPU-based general solvers with in-place LU decomposition + /// @tparam SolverParametersPolicy Parameters for the ODE solver + /// @tparam DenseMatrixPolicy Policy for dense matrices + /// @tparam SparseMatrixPolicy Policy for sparse matrices + /// @tparam LuDecompositionPolicy Policy for the LU decomposition + template< + class SolverParametersPolicy, + class DenseMatrix = Matrix, + class SparseMatrixPolicy = SparseMatrix, + class LuDecompositionPolicy = LuDecompositionInPlace> + using CpuSolverBuilderInPlace = SolverBuilder< + SolverParametersPolicy, + DenseMatrix, + SparseMatrixPolicy, + ProcessSet, + LuDecompositionPolicy, + LinearSolverInPlace, + State>; + } // namespace micm #include "solver_builder.inl" \ No newline at end of file diff --git a/test/integration/test_analytical_backward_euler.cpp b/test/integration/test_analytical_backward_euler.cpp index 7fb703806..5337c7213 100644 --- a/test/integration/test_analytical_backward_euler.cpp +++ b/test/integration/test_analytical_backward_euler.cpp @@ -15,8 +15,7 @@ using VectorBackwardEuler = micm::CpuSolverBuilder< micm::VectorMatrix, micm::SparseMatrix>>; template -using VectorStateType = - micm::State, micm::SparseMatrix>>; +using VectorStateType = typename VectorBackwardEuler::StatePolicyType; template using VectorBackwardEulerDoolittle = micm::CpuSolverBuilder< @@ -26,10 +25,7 @@ using VectorBackwardEulerDoolittle = micm::CpuSolverBuilder< micm::LuDecompositionDoolittle>; template -using VectorStateTypeDoolittle = micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionDoolittle>; +using VectorStateTypeDoolittle = typename VectorBackwardEulerDoolittle::StatePolicyType; template using VectorBackwardEulerDolittleCSC = micm::CpuSolverBuilder< @@ -39,10 +35,7 @@ using VectorBackwardEulerDolittleCSC = micm::CpuSolverBuilder< micm::LuDecompositionDoolittle>; template -using VectorStateTypeDoolittleCSC = micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionDoolittle>; +using VectorStateTypeDoolittleCSC = typename VectorBackwardEulerDolittleCSC::StatePolicyType; template using VectorBackwardEulerMozart = micm::CpuSolverBuilder< @@ -52,10 +45,7 @@ using VectorBackwardEulerMozart = micm::CpuSolverBuilder< micm::LuDecompositionMozart>; template -using VectorStateTypeMozart = micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionMozart>; +using VectorStateTypeMozart = typename VectorBackwardEulerMozart::StatePolicyType; template using VectorBackwardEulerMozartCSC = micm::CpuSolverBuilder< @@ -65,52 +55,27 @@ using VectorBackwardEulerMozartCSC = micm::CpuSolverBuilder< micm::LuDecompositionMozart>; template -using VectorStateTypeMozartCSC = micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionMozart>; +using VectorStateTypeMozartCSC = typename VectorBackwardEulerMozartCSC::StatePolicyType; template -using VectorBackwardEulerMozartInPlace = micm::SolverBuilder< +using VectorBackwardEulerDoolittleInPlace = micm::CpuSolverBuilderInPlace< micm::BackwardEulerSolverParameters, micm::VectorMatrix, micm::SparseMatrix>, - micm::ProcessSet, - micm::LuDecompositionMozartInPlace, - micm::LinearSolverInPlace< - micm::SparseMatrix>, - micm::LuDecompositionMozartInPlace>, - micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionMozartInPlace>>; + micm::LuDecompositionDoolittleInPlace>; template -using VectorStateTypeMozartInPlace = micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionMozartInPlace>; +using VectorStateTypeDoolittleInPlace = typename VectorBackwardEulerDoolittleInPlace::StatePolicyType; template -using VectorBackwardEulerMozartInPlace = micm::SolverBuilder< +using VectorBackwardEulerMozartInPlace = micm::CpuSolverBuilderInPlace< micm::BackwardEulerSolverParameters, micm::VectorMatrix, micm::SparseMatrix>, - micm::ProcessSet, - micm::LuDecompositionMozartInPlace, - micm::LinearSolverInPlace< - micm::SparseMatrix>, - micm::LuDecompositionMozartInPlace>, - micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionMozartInPlace>>; + micm::LuDecompositionMozartInPlace>; template -using VectorStateTypeMozartInPlace = micm::State< - micm::VectorMatrix, - micm::SparseMatrix>, - micm::LuDecompositionMozartInPlace>; +using VectorStateTypeMozartInPlace = typename VectorBackwardEulerMozartInPlace::StatePolicyType; auto backward_euler = micm::CpuSolverBuilder(micm::BackwardEulerSolverParameters()); auto backard_euler_vector_1 = VectorBackwardEuler<1>(micm::BackwardEulerSolverParameters()); @@ -134,11 +99,16 @@ auto backward_euler_vector_mozart_csc_1 = VectorBackwardEulerMozartCSC<1>(micm:: auto backward_euler_vector_mozart_csc_2 = VectorBackwardEulerMozartCSC<2>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_csc_3 = VectorBackwardEulerMozartCSC<3>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_csc_4 = VectorBackwardEulerMozartCSC<4>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_1 = VectorBackwardEulerDoolittleInPlace<1>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_2 = VectorBackwardEulerDoolittleInPlace<2>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_3 = VectorBackwardEulerDoolittleInPlace<3>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_4 = VectorBackwardEulerDoolittleInPlace<4>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_1 = VectorBackwardEulerMozartInPlace<1>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_2 = VectorBackwardEulerMozartInPlace<2>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_3 = VectorBackwardEulerMozartInPlace<3>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_4 = VectorBackwardEulerMozartInPlace<4>(micm::BackwardEulerSolverParameters()); + TEST(AnalyticalExamples, Troe) { test_analytical_troe(backward_euler, 1e-6); @@ -174,6 +144,14 @@ TEST(AnalyticalExamples, Troe) backward_euler_vector_mozart_csc_3, 1e-6); test_analytical_troe, VectorStateTypeMozartCSC<4>>( backward_euler_vector_mozart_csc_4, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleInPlace<1>>( + backward_euler_vector_doolittle_in_place_1, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleInPlace<2>>( + backward_euler_vector_doolittle_in_place_2, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleInPlace<3>>( + backward_euler_vector_doolittle_in_place_3, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleInPlace<4>>( + backward_euler_vector_doolittle_in_place_4, 1e-6); test_analytical_troe, VectorStateTypeMozartInPlace<1>>( backward_euler_vector_mozart_in_place_1, 1e-6); test_analytical_troe, VectorStateTypeMozartInPlace<2>>( diff --git a/test/unit/solver/CMakeLists.txt b/test/unit/solver/CMakeLists.txt index 826637172..b30e40a98 100644 --- a/test/unit/solver/CMakeLists.txt +++ b/test/unit/solver/CMakeLists.txt @@ -11,6 +11,7 @@ create_standard_test(NAME linear_solver SOURCES test_linear_solver.cpp) create_standard_test(NAME linear_solver_in_place SOURCES test_linear_solver_in_place.cpp) create_standard_test(NAME lu_decomposition_doolittle SOURCES test_lu_decomposition_doolittle.cpp) create_standard_test(NAME lu_decomposition_mozart SOURCES test_lu_decomposition_mozart.cpp) +create_standard_test(NAME lu_decomposition_doolittle_in_place SOURCES test_lu_decomposition_doolittle_in_place.cpp) create_standard_test(NAME lu_decomposition_mozart_in_place SOURCES test_lu_decomposition_mozart_in_place.cpp) create_standard_test(NAME rosenbrock SOURCES test_rosenbrock.cpp) create_standard_test(NAME state SOURCES test_state.cpp) diff --git a/test/unit/solver/test_lu_decomposition_doolittle_in_place.cpp b/test/unit/solver/test_lu_decomposition_doolittle_in_place.cpp new file mode 100644 index 000000000..3156ae080 --- /dev/null +++ b/test/unit/solver/test_lu_decomposition_doolittle_in_place.cpp @@ -0,0 +1,75 @@ +#include "test_lu_decomposition_in_place_policy.hpp" + +#include +#include +#include + +#include + +using SparseMatrixTest = micm::SparseMatrix; + +using Group1SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group3SparseVectorMatrix = micm::SparseMatrix>; +using Group4SparseVectorMatrix = micm::SparseMatrix>; + +TEST(LuDecompositionDoolittleInPlace, DenseMatrixStandardOrdering) +{ + testDenseMatrix(); +} + +TEST(LuDecompositionDoolittleInPlace, RandomMatrixStandardOrdering) +{ + testRandomMatrix(1); + testRandomMatrix(5); +} + +TEST(LuDecompositionDoolittleInPlace, DiagonalMatrixStandardOrdering) +{ + testDiagonalMatrix(5); +} + +TEST(LuDecompositionDoolittleInPlace, AgnosticToInitialValueStandardOrdering) +{ + double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; + for (auto& value : initial_values) + { + testExtremeValueInitialization(5, value); + } +} + +TEST(LuDecompositionDoolittleInPlace, DenseMatrixVectorOrdering) +{ + testDenseMatrix(); + testDenseMatrix(); + testDenseMatrix(); + testDenseMatrix(); +} + +TEST(LuDecompositionDoolittleInPlace, RandomMatrixVectorOrdering) +{ + testRandomMatrix(5); + testRandomMatrix(5); + testRandomMatrix(5); + testRandomMatrix(5); +} + +TEST(LuDecompositionDoolittleInPlace, DiagonalMatrixVectorOrdering) +{ + testDiagonalMatrix(5); + testDiagonalMatrix(5); + testDiagonalMatrix(5); + testDiagonalMatrix(5); +} + +TEST(LuDecompositionDoolittleInPlace, AgnosticToInitialValueVectorOrdering) +{ + double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; + for (auto& value : initial_values) + { + testExtremeValueInitialization(5, value); + testExtremeValueInitialization(5, value); + testExtremeValueInitialization(5, value); + testExtremeValueInitialization(5, value); + } +} \ No newline at end of file From faf21bd72c8dd8b57b780dff55016f54177cd1d3 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 9 Jan 2025 15:20:23 -0800 Subject: [PATCH 2/2] Auto-format code changes (#706) Auto-format code using Clang-Format Co-authored-by: GitHub Actions --- include/micm/solver/lu_decomposition.hpp | 2 +- .../solver/lu_decomposition_doolittle_in_place.inl | 8 ++++---- include/micm/solver/solver_builder.hpp | 2 +- test/integration/test_analytical_backward_euler.cpp | 13 ++++++++----- 4 files changed, 14 insertions(+), 11 deletions(-) diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 9c04fb8c8..66861fea7 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -3,8 +3,8 @@ #pragma once #include "lu_decomposition_doolittle.hpp" -#include "lu_decomposition_mozart.hpp" #include "lu_decomposition_doolittle_in_place.hpp" +#include "lu_decomposition_mozart.hpp" #include "lu_decomposition_mozart_in_place.hpp" #include diff --git a/include/micm/solver/lu_decomposition_doolittle_in_place.inl b/include/micm/solver/lu_decomposition_doolittle_in_place.inl index 514ead06c..93535029c 100644 --- a/include/micm/solver/lu_decomposition_doolittle_in_place.inl +++ b/include/micm/solver/lu_decomposition_doolittle_in_place.inl @@ -179,16 +179,16 @@ namespace micm const std::size_t ALU_BlockSize = ALU.NumberOfBlocks(); constexpr std::size_t ALU_GroupVectorSize = SparseMatrixPolicy::GroupVectorSize(); const std::size_t ALU_GroupSizeOfFlatBlockSize = ALU.GroupSize(); - + // Loop over groups of blocks for (std::size_t i_group = 0; i_group < ALU.NumberOfGroups(ALU_BlockSize); ++i_group) { - auto ALU_vector = std::next(ALU.AsVector().begin(), i_group * ALU_GroupSizeOfFlatBlockSize); + auto ALU_vector = std::next(ALU.AsVector().begin(), i_group * ALU_GroupSizeOfFlatBlockSize); const std::size_t n_cells = std::min(ALU_GroupVectorSize, ALU_BlockSize - i_group * ALU_GroupVectorSize); auto aik_njk = aik_njk_.begin(); auto aij_ajk = aij_ajk_.begin(); auto aki_nji = aki_nji_.begin(); - auto akj_aji = akj_aji_.begin(); + auto akj_aji = akj_aji_.begin(); for (const auto& nik_nki_aii : nik_nki_aii_) { for (std::size_t ik = 0; ik < std::get<0>(nik_nki_aii); ++ik) @@ -213,7 +213,7 @@ namespace micm ALU_vector[aki_nji->first + i] /= ALU_vector[std::get<2>(nik_nki_aii) + i]; ++aki_nji; } - } + } } } diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index 0f6c33f2f..066871bc8 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -137,7 +137,7 @@ namespace micm State>; /// @brief Builder of CPU-based general solvers with in-place LU decomposition - /// @tparam SolverParametersPolicy Parameters for the ODE solver + /// @tparam SolverParametersPolicy Parameters for the ODE solver /// @tparam DenseMatrixPolicy Policy for dense matrices /// @tparam SparseMatrixPolicy Policy for sparse matrices /// @tparam LuDecompositionPolicy Policy for the LU decomposition diff --git a/test/integration/test_analytical_backward_euler.cpp b/test/integration/test_analytical_backward_euler.cpp index 5337c7213..9ed9653ad 100644 --- a/test/integration/test_analytical_backward_euler.cpp +++ b/test/integration/test_analytical_backward_euler.cpp @@ -99,16 +99,19 @@ auto backward_euler_vector_mozart_csc_1 = VectorBackwardEulerMozartCSC<1>(micm:: auto backward_euler_vector_mozart_csc_2 = VectorBackwardEulerMozartCSC<2>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_csc_3 = VectorBackwardEulerMozartCSC<3>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_csc_4 = VectorBackwardEulerMozartCSC<4>(micm::BackwardEulerSolverParameters()); -auto backward_euler_vector_doolittle_in_place_1 = VectorBackwardEulerDoolittleInPlace<1>(micm::BackwardEulerSolverParameters()); -auto backward_euler_vector_doolittle_in_place_2 = VectorBackwardEulerDoolittleInPlace<2>(micm::BackwardEulerSolverParameters()); -auto backward_euler_vector_doolittle_in_place_3 = VectorBackwardEulerDoolittleInPlace<3>(micm::BackwardEulerSolverParameters()); -auto backward_euler_vector_doolittle_in_place_4 = VectorBackwardEulerDoolittleInPlace<4>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_1 = + VectorBackwardEulerDoolittleInPlace<1>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_2 = + VectorBackwardEulerDoolittleInPlace<2>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_3 = + VectorBackwardEulerDoolittleInPlace<3>(micm::BackwardEulerSolverParameters()); +auto backward_euler_vector_doolittle_in_place_4 = + VectorBackwardEulerDoolittleInPlace<4>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_1 = VectorBackwardEulerMozartInPlace<1>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_2 = VectorBackwardEulerMozartInPlace<2>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_3 = VectorBackwardEulerMozartInPlace<3>(micm::BackwardEulerSolverParameters()); auto backward_euler_vector_mozart_in_place_4 = VectorBackwardEulerMozartInPlace<4>(micm::BackwardEulerSolverParameters()); - TEST(AnalyticalExamples, Troe) { test_analytical_troe(backward_euler, 1e-6);