From 7eccbb8f7b1a918b6d432dceedcf7b0f37906c57 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 10:46:02 -0700 Subject: [PATCH 01/21] Add dense LA solvers. --- .../denseLinearAlgebra/CMakeLists.txt | 1 + .../denseLinearAlgebra/denseLASolvers.hpp | 255 ++++++++++++++++++ 2 files changed, 256 insertions(+) create mode 100644 src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp diff --git a/src/coreComponents/denseLinearAlgebra/CMakeLists.txt b/src/coreComponents/denseLinearAlgebra/CMakeLists.txt index 791767bc5c6..55380d8c4a5 100644 --- a/src/coreComponents/denseLinearAlgebra/CMakeLists.txt +++ b/src/coreComponents/denseLinearAlgebra/CMakeLists.txt @@ -1,6 +1,7 @@ # Specify all headers set( denseLinearAlgebra_headers common/layouts.hpp + denseLASolvers.hpp interfaces/blaslapack/BlasLapackFunctions.h interfaces/blaslapack/BlasLapackLA.hpp ) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp new file mode 100644 index 00000000000..793820cd779 --- /dev/null +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -0,0 +1,255 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file denseLASolvers.hpp + */ +#ifndef GEOS_DENSELINEARALGEBRA_DENSELASOLVERS_HPP_ +#define GEOS_DENSELINEARALGEBRA_DENSELASOLVERS_HPP_ + +#include "common/DataTypes.hpp" +#include "denseLinearAlgebra/common/layouts.hpp" + +#include + +namespace geos +{ + +namespace denseLinearAlgebra +{ + +namespace internal +{ +/** + * @brief Solves a 2x2 linear system A * x = b. + * + * This function solves a linear system of the form A * x = b, where A is a 2x2 matrix, + * b is a 2x1 vector, and x is the solution vector. The function checks the sizes + * of the inputs to ensure they conform to the expected dimensions. It also checks that + * the determinant of matrix A is not near zero to avoid solving a singular system. + * + * @tparam MATRIX_TYPE The type of the matrix A. Must support indexing with `A[i][j]`. + * @tparam RHS_TYPE The type of the right-hand side vector b. Must support indexing with `b[i]`. + * @tparam SOL_TYPE The type of the solution vector x. Must support indexing with `x[i]`. + * + * @param[in] A The 2x2 matrix representing the system of equations. Must have size 2x2. + * @param[in] b The 2-element vector representing the right-hand side of the equation. + * @param[out] x The 2-element vector that will store the solution to the system. + */ +template< typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +GEOS_HOST_DEVICE +inline +void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x) +{ + LvArray::tensorOps::internal::checkSizes< 2, 2 >( A ); + LvArray::tensorOps::internal::checkSizes< 2 >( b ); + LvArray::tensorOps::internal::checkSizes< 2 >( x ); + + real64 const detA = A[0][0] * A[1][1] - A[0][1] * A[1][0]; + + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs(detA), LvArray::NumericLimits< real64 >::epsilon;, "Singular system." ); + + x[0] = (A[1][1] * b[0] - A[0][1] * b[1] ) / detA; + x[1] = (A[0][0] * b[1] - A[1][0] * b[0] ) / detA; +} + +/** + * @brief Solves a 3x3 linear system A * x = b. + * + * This function solves a linear system of the form A * x = b, where A is a 3x3 matrix, + * b is a 3x1 vector, and x is the solution vector. The function checks the sizes + * of the inputs to ensure they conform to the expected dimensions. It also checks that + * the determinant of matrix A is not near zero to avoid solving a singular system. + * + * @tparam MATRIX_TYPE The type of the matrix A. Must support indexing with `A[i][j]`. + * @tparam RHS_TYPE The type of the right-hand side vector b. Must support indexing with `b[i]`. + * @tparam SOL_TYPE The type of the solution vector x. Must support indexing with `x[i]`. + * + * @param[in] A The 3x3 matrix representing the system of equations. Must have size 3x3. + * @param[in] b The 3-element vector representing the right-hand side of the equation. + * @param[out] x The 3-element vector that will store the solution to the system. + */ +template< typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +GEOS_HOST_DEVICE +inline +void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x) +{ + LvArray::tensorOps::internal::checkSizes< 3, 3 >( A ); + LvArray::tensorOps::internal::checkSizes< 3 >( b ); + LvArray::tensorOps::internal::checkSizes< 3 >( x ); + + real64 const detA = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - + A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) + + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]); + + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs(detA), LvArray::NumericLimits< real64 >::epsilon;, "Singular system." ); + + real64 const detX0 = b[0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - + b[1] * (A[0][1] * A[2][2] - A[0][2] * A[2][1]) + + b[2] * (A[0][1] * A[1][2] - A[0][2] * A[1][1]); + + real64 const detX1 = A[0][0] * (b[1] * A[2][2] - b[2] * A[2][1]) - + A[0][1] * (b[0] * A[2][2] - b[2] * A[2][0]) + + A[0][2] * (b[0] * A[1][2] - b[1] * A[1][0]); + + real64 const detX2 = A[0][0] * (A[1][1] * b[2] - A[1][2] * b[1]) - + A[0][1] * (A[1][0] * b[2] - A[1][2] * b[0]) + + A[0][2] * (A[1][0] * b[1] - A[1][1] * b[0]); + + x[0] = detX0 / detA; + x[1] = detX1 / detA; + x[2] = detX2 / detA; +} + + +template< std::ptrdiff_t N, + typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE & b, SOL_TYPE && x ) +{ + for (std::ptrdiff_t i = N - 1; i >= 0; --i) + { + real64 sum = b[i]; + for (std::ptrdiff_t j = i + 1; j < N; ++j) + { + sum -= A[i][j] * x[j]; + } + x[i] = sum / A[i][i]; + } +} + +template< std::ptrdiff_t N, + typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +GEOS_HOST_DEVICE +inline +void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +{ + static_assert( N > 0, "N must be greater than 0." ); + internal::checkSizes< N, N >( matrix ); + internal::checkSizes< N >( b ); + internal::checkSizes< N >( x ); + + + // Step 1: Transform into an upper triangular matrix + + // 1.a. Find the pivot + for (std::ptrdiff_t i = 0; i < N; ++i) + { + std::ptrdiff_t max_row = i; + for (std::ptrdiff_t k = i + 1; k < N; ++k) + { + if (std::abs(A[k][i]) > std::abs(A[max_row][i])) + { + max_row = k; + } + } + + // 1.b. Swap rows + for (std::ptrdiff_t k = i; k < N; ++k) + { + std::swap(A[i][k], A[max_row][k]); + } + std::swap(b[i], b[max_row]); + + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs(A[i][i]), LvArray::NumericLimits< real64 >::epsilon, "Singular matrix." ); + + // 1.c Eliminate entries below the pivot + for (std::ptrdiff_t k = i + 1; k < N; ++k) + { + real64 const scaling = A[k][i] / A[i][i]; + for (std::ptrdiff_t j = i; j < N; ++j) + { + A[k][j] -= scaling * A[i][j]; + } + b[k] -= scaling * b[i]; + } + } + + // Step 2: Backward substitution + solveUpperTriangularSystem( A, b, std::forward(x) ) +} + +/** + * Const version of the function + */ +template< std::ptrdiff_t N, + typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +GEOS_HOST_DEVICE +inline +void solveGaussianElimination( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +{ + real64[N][N] A_copy{}; + real64[N] b_copy{}; + + for(std::ptrdiff_t i=0; i < N; ++j) + { + b_copy[i] = b[i]; + for( std::ptrdiff_t j=0; j < N; ++j ) + { + A_copy[i][j] = A[i][j]; + }; + }; + + solveGaussianElimination( A_copy, b_copy, std::forward(x) ); +} + +}; // internal namespace + + +/** + * + */ +template< std::ptrdiff_t N, + typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +GEOS_HOST_DEVICE +inline +void solve( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +{ + static_assert( N > 0, "N must be greater than 0." ); + internal::checkSizes< N, N >( A ); + internal::checkSizes< N >( b ); + internal::checkSizes< N >( x ); + + if constexpr ( N == 2 ) + { + internal::solveTwoByTwoSystem( A, b, std::forward(x) ); + } + else if constexpr( N == 3 ) + { + internal::solveThreeByThreeSystem( A, b, std::forward(x) ); + } + else + { + internal::solveGaussianElimination< N >( A, b, std::forward(x) ); + } +} + +}; + +}; + + +#endif /*GEOS_DENSELINEARALGEBRA_DENSELASOLVERS_HPP_*/ From d3d6b83fea28d0678edfd5350a0d8813bdf3cc7c Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 11:00:08 -0700 Subject: [PATCH 02/21] add documentation. --- .../denseLinearAlgebra/denseLASolvers.hpp | 56 ++++++++++++++++++- 1 file changed, 53 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 793820cd779..5955730fe04 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -117,7 +117,20 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP x[2] = detX2 / detA; } - +/** + * @brief Solves a linear system where the matrix is upper triangular using back substitution. + * + * This function solves the linear system `Ax = b`, where `A` is an upper triangular matrix, using + * back substitution. The solution `x` is computed and stored in the provided output vector. + * + * @tparam N The size of the square matrix `A`. + * @tparam MATRIX_TYPE The type of the matrix `A`. + * @tparam RHS_TYPE The type of the right-hand side vector `b`. + * @tparam SOL_TYPE The type of the solution vector `x`. + * @param[in] A The upper triangular matrix representing the coefficients of the system. + * @param[in,out] b The right-hand side vector. It is used to compute the solution, but the values are not preserved. + * @param[out] x The solution vector. The result of solving the system `Ax = b` using back substitution. + */ template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, @@ -135,6 +148,21 @@ void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE & b, SOL_TYPE & } } +/** + * @brief Solves a linear system using Gaussian elimination. + * + * This function performs Gaussian elimination on the given matrix `A` and right-hand side vector `b`. + * It transforms the matrix `A` into an upper triangular matrix and then solves for the solution `x` + * using back substitution. + * + * @tparam N The size of the square matrix `A`. + * @tparam MATRIX_TYPE The type of the matrix `A`. + * @tparam RHS_TYPE The type of the right-hand side vector `b`. + * @tparam SOL_TYPE The type of the solution vector `x`. + * @param[in,out] A The matrix to be transformed into an upper triangular matrix. Modified in place. + * @param[in,out] b The right-hand side vector. Modified in place to reflect the transformed system. + * @param[out] x The solution vector. The result of solving the system `Ax = b`. + */ template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, @@ -189,7 +217,18 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) } /** - * Const version of the function + * @brief Const version of the solveGaussianElimination function. + * + * This function solves a linear system using Gaussian elimination, without modifying the original matrix `A` + * and vector `b`. It creates copies of `A` and `b`, performs the elimination, and then solves the system. + * + * @tparam N The size of the square matrix `A`. + * @tparam MATRIX_TYPE The type of the matrix `A`. + * @tparam RHS_TYPE The type of the right-hand side vector `b`. + * @tparam SOL_TYPE The type of the solution vector `x`. + * @param[in] A The constant matrix representing the coefficients of the system. + * @param[in] b The constant right-hand side vector. + * @param[out] x The solution vector. The result of solving the system `Ax = b`. */ template< std::ptrdiff_t N, typename MATRIX_TYPE, @@ -216,9 +255,20 @@ void solveGaussianElimination( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TY }; // internal namespace - /** + * @brief Solves a linear system using the most appropriate method based on the size of the system. + * + * This function determines the appropriate method for solving a linear system `Ax = b` based on + * the size of the matrix `A`. For 2x2 and 3x3 systems, specialized solvers are used. For larger systems, + * Gaussian elimination is employed. * + * @tparam N The size of the square matrix `A`. + * @tparam MATRIX_TYPE The type of the matrix `A`. + * @tparam RHS_TYPE The type of the right-hand side vector `b`. + * @tparam SOL_TYPE The type of the solution vector `x`. + * @param[in] A The constant matrix representing the coefficients of the system. + * @param[in] b The constant right-hand side vector. + * @param[out] x The solution vector. The result of solving the system `Ax = b`. */ template< std::ptrdiff_t N, typename MATRIX_TYPE, From d3c693e2f43b1d7729a3fe82506d53ae236dcccd Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 11:07:28 -0700 Subject: [PATCH 03/21] add non-const interface. --- .../denseLinearAlgebra/denseLASolvers.hpp | 46 ++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 5955730fe04..c6bbfa1c6bb 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -260,7 +260,50 @@ void solveGaussianElimination( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TY * * This function determines the appropriate method for solving a linear system `Ax = b` based on * the size of the matrix `A`. For 2x2 and 3x3 systems, specialized solvers are used. For larger systems, - * Gaussian elimination is employed. + * Gaussian elimination is employed. The matrix and the rhs are modified by the function. + * + * @tparam N The size of the square matrix `A`. + * @tparam MATRIX_TYPE The type of the matrix `A`. + * @tparam RHS_TYPE The type of the right-hand side vector `b`. + * @tparam SOL_TYPE The type of the solution vector `x`. + * @param[in] A The constant matrix representing the coefficients of the system. + * @param[in] b The constant right-hand side vector. + * @param[out] x The solution vector. The result of solving the system `Ax = b`. + */ +template< std::ptrdiff_t N, + typename MATRIX_TYPE, + typename RHS_TYPE, + typename SOL_TYPE > +GEOS_HOST_DEVICE +inline +void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +{ + static_assert( N > 0, "N must be greater than 0." ); + internal::checkSizes< N, N >( A ); + internal::checkSizes< N >( b ); + internal::checkSizes< N >( x ); + + if constexpr ( N == 2 ) + { + internal::solveTwoByTwoSystem( A, b, std::forward(x) ); + } + else if constexpr( N == 3 ) + { + internal::solveThreeByThreeSystem( A, b, std::forward(x) ); + } + else + { + internal::solveGaussianElimination< N >( A, b, std::forward(x) ); + } +} + +/** + * @brief Solves a linear system using the most appropriate method based on the size of the system. + * + * This function determines the appropriate method for solving a linear system `Ax = b` based on + * the size of the matrix `A`. For 2x2 and 3x3 systems, specialized solvers are used. For larger systems, + * Gaussian elimination is employed. The matrix `A` and the right-hand side vector `b` are not modified, + * making this function suitable for cases where `A` and `b` need to remain constant. * * @tparam N The size of the square matrix `A`. * @tparam MATRIX_TYPE The type of the matrix `A`. @@ -297,6 +340,7 @@ void solve( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) } } + }; }; From f855e3fbdf9a7c29d2ec48943e44bdbf5faa89de Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 11:53:33 -0700 Subject: [PATCH 04/21] suggestions from code review --- .../denseLinearAlgebra/denseLASolvers.hpp | 251 +++++++----------- 1 file changed, 95 insertions(+), 156 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index c6bbfa1c6bb..bf714c1303b 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -31,13 +31,13 @@ namespace denseLinearAlgebra { namespace internal -{ +{ /** * @brief Solves a 2x2 linear system A * x = b. * - * This function solves a linear system of the form A * x = b, where A is a 2x2 matrix, - * b is a 2x1 vector, and x is the solution vector. The function checks the sizes - * of the inputs to ensure they conform to the expected dimensions. It also checks that + * This function solves a linear system of the form A * x = b, where A is a 2x2 matrix, + * b is a 2x1 vector, and x is the solution vector. The function checks the sizes + * of the inputs to ensure they conform to the expected dimensions. It also checks that * the determinant of matrix A is not near zero to avoid solving a singular system. * * @tparam MATRIX_TYPE The type of the matrix A. Must support indexing with `A[i][j]`. @@ -48,31 +48,31 @@ namespace internal * @param[in] b The 2-element vector representing the right-hand side of the equation. * @param[out] x The 2-element vector that will store the solution to the system. */ -template< typename MATRIX_TYPE, - typename RHS_TYPE, +template< typename MATRIX_TYPE, + typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x) +void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { - LvArray::tensorOps::internal::checkSizes< 2, 2 >( A ); - LvArray::tensorOps::internal::checkSizes< 2 >( b ); - LvArray::tensorOps::internal::checkSizes< 2 >( x ); + LvArray::tensorOps::internal::checkSizes< 2, 2 >( A ); + LvArray::tensorOps::internal::checkSizes< 2 >( b ); + LvArray::tensorOps::internal::checkSizes< 2 >( x ); - real64 const detA = A[0][0] * A[1][1] - A[0][1] * A[1][0]; + real64 const detA = A[0][0] * A[1][1] - A[0][1] * A[1][0]; - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs(detA), LvArray::NumericLimits< real64 >::epsilon;, "Singular system." ); + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon; , "Singular system." ); - x[0] = (A[1][1] * b[0] - A[0][1] * b[1] ) / detA; - x[1] = (A[0][0] * b[1] - A[1][0] * b[0] ) / detA; + x[0] = (A[1][1] * b[0] - A[0][1] * b[1] ) / detA; + x[1] = (A[0][0] * b[1] - A[1][0] * b[0] ) / detA; } /** * @brief Solves a 3x3 linear system A * x = b. * - * This function solves a linear system of the form A * x = b, where A is a 3x3 matrix, - * b is a 3x1 vector, and x is the solution vector. The function checks the sizes - * of the inputs to ensure they conform to the expected dimensions. It also checks that + * This function solves a linear system of the form A * x = b, where A is a 3x3 matrix, + * b is a 3x1 vector, and x is the solution vector. The function checks the sizes + * of the inputs to ensure they conform to the expected dimensions. It also checks that * the determinant of matrix A is not near zero to avoid solving a singular system. * * @tparam MATRIX_TYPE The type of the matrix A. Must support indexing with `A[i][j]`. @@ -83,46 +83,46 @@ void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && * @param[in] b The 3-element vector representing the right-hand side of the equation. * @param[out] x The 3-element vector that will store the solution to the system. */ -template< typename MATRIX_TYPE, - typename RHS_TYPE, +template< typename MATRIX_TYPE, + typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x) +void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { - LvArray::tensorOps::internal::checkSizes< 3, 3 >( A ); - LvArray::tensorOps::internal::checkSizes< 3 >( b ); - LvArray::tensorOps::internal::checkSizes< 3 >( x ); + LvArray::tensorOps::internal::checkSizes< 3, 3 >( A ); + LvArray::tensorOps::internal::checkSizes< 3 >( b ); + LvArray::tensorOps::internal::checkSizes< 3 >( x ); - real64 const detA = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) + - A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]); + real64 const detA = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - + A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) + + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs(detA), LvArray::NumericLimits< real64 >::epsilon;, "Singular system." ); + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon; , "Singular system." ); - real64 const detX0 = b[0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - - b[1] * (A[0][1] * A[2][2] - A[0][2] * A[2][1]) + - b[2] * (A[0][1] * A[1][2] - A[0][2] * A[1][1]); + real64 const detX0 = b[0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - + b[1] * (A[0][1] * A[2][2] - A[0][2] * A[2][1]) + + b[2] * (A[0][1] * A[1][2] - A[0][2] * A[1][1]); - real64 const detX1 = A[0][0] * (b[1] * A[2][2] - b[2] * A[2][1]) - - A[0][1] * (b[0] * A[2][2] - b[2] * A[2][0]) + - A[0][2] * (b[0] * A[1][2] - b[1] * A[1][0]); + real64 const detX1 = A[0][0] * (b[1] * A[2][2] - b[2] * A[2][1]) - + A[0][1] * (b[0] * A[2][2] - b[2] * A[2][0]) + + A[0][2] * (b[0] * A[1][2] - b[1] * A[1][0]); - real64 const detX2 = A[0][0] * (A[1][1] * b[2] - A[1][2] * b[1]) - - A[0][1] * (A[1][0] * b[2] - A[1][2] * b[0]) + - A[0][2] * (A[1][0] * b[1] - A[1][1] * b[0]); + real64 const detX2 = A[0][0] * (A[1][1] * b[2] - A[1][2] * b[1]) - + A[0][1] * (A[1][0] * b[2] - A[1][2] * b[0]) + + A[0][2] * (A[1][0] * b[1] - A[1][1] * b[0]); - x[0] = detX0 / detA; - x[1] = detX1 / detA; - x[2] = detX2 / detA; + x[0] = detX0 / detA; + x[1] = detX1 / detA; + x[2] = detX2 / detA; } /** * @brief Solves a linear system where the matrix is upper triangular using back substitution. - * - * This function solves the linear system `Ax = b`, where `A` is an upper triangular matrix, using + * + * This function solves the linear system `Ax = b`, where `A` is an upper triangular matrix, using * back substitution. The solution `x` is computed and stored in the provided output vector. - * + * * @tparam N The size of the square matrix `A`. * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. @@ -131,16 +131,16 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP * @param[in,out] b The right-hand side vector. It is used to compute the solution, but the values are not preserved. * @param[out] x The solution vector. The result of solving the system `Ax = b` using back substitution. */ -template< std::ptrdiff_t N, +template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE & b, SOL_TYPE && x ) { - for (std::ptrdiff_t i = N - 1; i >= 0; --i) + for( std::ptrdiff_t i = N - 1; i >= 0; --i ) { real64 sum = b[i]; - for (std::ptrdiff_t j = i + 1; j < N; ++j) + for( std::ptrdiff_t j = i + 1; j < N; ++j ) { sum -= A[i][j] * x[j]; } @@ -150,11 +150,11 @@ void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE & b, SOL_TYPE & /** * @brief Solves a linear system using Gaussian elimination. - * + * * This function performs Gaussian elimination on the given matrix `A` and right-hand side vector `b`. * It transforms the matrix `A` into an upper triangular matrix and then solves for the solution `x` * using back substitution. - * + * * @tparam N The size of the square matrix `A`. * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. @@ -163,48 +163,48 @@ void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE & b, SOL_TYPE & * @param[in,out] b The right-hand side vector. Modified in place to reflect the transformed system. * @param[out] x The solution vector. The result of solving the system `Ax = b`. */ -template< std::ptrdiff_t N, +template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > -GEOS_HOST_DEVICE +GEOS_HOST_DEVICE inline void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); internal::checkSizes< N, N >( matrix ); - internal::checkSizes< N >( b ); + internal::checkSizes< N >( b ); internal::checkSizes< N >( x ); - - // Step 1: Transform into an upper triangular matrix - + + // Step 1: Transform into an upper triangular matrix + // 1.a. Find the pivot - for (std::ptrdiff_t i = 0; i < N; ++i) + for( std::ptrdiff_t i = 0; i < N; ++i ) { std::ptrdiff_t max_row = i; - for (std::ptrdiff_t k = i + 1; k < N; ++k) + for( std::ptrdiff_t k = i + 1; k < N; ++k ) { - if (std::abs(A[k][i]) > std::abs(A[max_row][i])) + if( std::abs( A[k][i] ) > std::abs( A[max_row][i] )) { max_row = k; } } // 1.b. Swap rows - for (std::ptrdiff_t k = i; k < N; ++k) + for( std::ptrdiff_t k = i; k < N; ++k ) { - std::swap(A[i][k], A[max_row][k]); + std::swap( A[i][k], A[max_row][k] ); } - std::swap(b[i], b[max_row]); + std::swap( b[i], b[max_row] ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs(A[i][i]), LvArray::NumericLimits< real64 >::epsilon, "Singular matrix." ); + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( A[i][i] ), LvArray::NumericLimits< real64 >::epsilon, "Singular matrix." ); // 1.c Eliminate entries below the pivot - for (std::ptrdiff_t k = i + 1; k < N; ++k) + for( std::ptrdiff_t k = i + 1; k < N; ++k ) { real64 const scaling = A[k][i] / A[i][i]; - for (std::ptrdiff_t j = i; j < N; ++j) + for( std::ptrdiff_t j = i; j < N; ++j ) { A[k][j] -= scaling * A[i][j]; } @@ -213,55 +213,18 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) } // Step 2: Backward substitution - solveUpperTriangularSystem( A, b, std::forward(x) ) -} - -/** - * @brief Const version of the solveGaussianElimination function. - * - * This function solves a linear system using Gaussian elimination, without modifying the original matrix `A` - * and vector `b`. It creates copies of `A` and `b`, performs the elimination, and then solves the system. - * - * @tparam N The size of the square matrix `A`. - * @tparam MATRIX_TYPE The type of the matrix `A`. - * @tparam RHS_TYPE The type of the right-hand side vector `b`. - * @tparam SOL_TYPE The type of the solution vector `x`. - * @param[in] A The constant matrix representing the coefficients of the system. - * @param[in] b The constant right-hand side vector. - * @param[out] x The solution vector. The result of solving the system `Ax = b`. - */ -template< std::ptrdiff_t N, - typename MATRIX_TYPE, - typename RHS_TYPE, - typename SOL_TYPE > -GEOS_HOST_DEVICE -inline -void solveGaussianElimination( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) -{ - real64[N][N] A_copy{}; - real64[N] b_copy{}; - - for(std::ptrdiff_t i=0; i < N; ++j) - { - b_copy[i] = b[i]; - for( std::ptrdiff_t j=0; j < N; ++j ) - { - A_copy[i][j] = A[i][j]; - }; - }; - - solveGaussianElimination( A_copy, b_copy, std::forward(x) ); + solveUpperTriangularSystem< N >( A, b, std::forward< N >( x ) ) } }; // internal namespace /** * @brief Solves a linear system using the most appropriate method based on the size of the system. - * - * This function determines the appropriate method for solving a linear system `Ax = b` based on - * the size of the matrix `A`. For 2x2 and 3x3 systems, specialized solvers are used. For larger systems, + * + * This function determines the appropriate method for solving a linear system `Ax = b` based on + * the size of the matrix `A`. For 2x2 and 3x3 systems, specialized solvers are used. For larger systems, * Gaussian elimination is employed. The matrix and the rhs are modified by the function. - * + * * @tparam N The size of the square matrix `A`. * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. @@ -270,77 +233,53 @@ void solveGaussianElimination( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TY * @param[in] b The constant right-hand side vector. * @param[out] x The solution vector. The result of solving the system `Ax = b`. */ -template< std::ptrdiff_t N, +template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, - typename SOL_TYPE > -GEOS_HOST_DEVICE + typename SOL_TYPE, + bool MODIFY_MATRIX = true > +GEOS_HOST_DEVICE inline void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); + static_assert( N < 10, "N cannot be larger than 9" ); internal::checkSizes< N, N >( A ); - internal::checkSizes< N >( b ); + internal::checkSizes< N >( b ); internal::checkSizes< N >( x ); - - if constexpr ( N == 2 ) - { - internal::solveTwoByTwoSystem( A, b, std::forward(x) ); - } - else if constexpr( N == 3 ) - { - internal::solveThreeByThreeSystem( A, b, std::forward(x) ); - } - else - { - internal::solveGaussianElimination< N >( A, b, std::forward(x) ); - } -} -/** - * @brief Solves a linear system using the most appropriate method based on the size of the system. - * - * This function determines the appropriate method for solving a linear system `Ax = b` based on - * the size of the matrix `A`. For 2x2 and 3x3 systems, specialized solvers are used. For larger systems, - * Gaussian elimination is employed. The matrix `A` and the right-hand side vector `b` are not modified, - * making this function suitable for cases where `A` and `b` need to remain constant. - * - * @tparam N The size of the square matrix `A`. - * @tparam MATRIX_TYPE The type of the matrix `A`. - * @tparam RHS_TYPE The type of the right-hand side vector `b`. - * @tparam SOL_TYPE The type of the solution vector `x`. - * @param[in] A The constant matrix representing the coefficients of the system. - * @param[in] b The constant right-hand side vector. - * @param[out] x The solution vector. The result of solving the system `Ax = b`. - */ -template< std::ptrdiff_t N, - typename MATRIX_TYPE, - typename RHS_TYPE, - typename SOL_TYPE > -GEOS_HOST_DEVICE -inline -void solve( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) -{ - static_assert( N > 0, "N must be greater than 0." ); - internal::checkSizes< N, N >( A ); - internal::checkSizes< N >( b ); - internal::checkSizes< N >( x ); - if constexpr ( N == 2 ) { - internal::solveTwoByTwoSystem( A, b, std::forward(x) ); + internal::solveTwoByTwoSystem( A, b, std::forward< SOL_TYPE >( x ) ); } - else if constexpr( N == 3 ) + else if constexpr ( N == 3 ) { - internal::solveThreeByThreeSystem( A, b, std::forward(x) ); + internal::solveThreeByThreeSystem( A, b, std::forward< SOL_TYPE >( x ) ); } else { - internal::solveGaussianElimination< N >( A, b, std::forward(x) ); + if constexpr ( MODIFY_MATRIX ) + { + internal::solveGaussianElimination< N >( A, b, std::forward< SOL_TYPE >( x ) ); + } + else + { + real64[N][N] A_copy{}; + real64[N] b_copy{}; + + for( std::ptrdiff_t i=0; i < N; ++j ) + { + b_copy[i] = b[i]; + for( std::ptrdiff_t j=0; j < N; ++j ) + { + A_copy[i][j] = A[i][j]; + } + } + internal::solveGaussianElimination< N >( A_copy, b_copy, std::forward< SOL_TYPE >( x ) ); + } } } - }; }; From 026017a39c5498a09590e04a01fef78e5107a078 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 11:56:45 -0700 Subject: [PATCH 05/21] adjust comments. --- src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index bf714c1303b..a727dea4567 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -229,6 +229,9 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. * @tparam SOL_TYPE The type of the solution vector `x`. + * @tparam MODIFY_MATRIX Boolean flag indicating whether the input matrix `A` and vector `b` should be modified. + * If `true`, the matrix `A` and vector `b` are modified in place. If `false`, copies of + * `A` and `b` are made, and the original data is left unchanged. * @param[in] A The constant matrix representing the coefficients of the system. * @param[in] b The constant right-hand side vector. * @param[out] x The solution vector. The result of solving the system `Ax = b`. From 640d56b46de90cc5c5c5ac6c671c7ee2fcf05922 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 21:00:07 -0700 Subject: [PATCH 06/21] fix bugs and add first two unitTests. --- .../denseLinearAlgebra/denseLASolvers.hpp | 82 ++++++------- .../unitTests/CMakeLists.txt | 3 +- .../unitTests/testDenseLASolvers.cpp | 108 ++++++++++++++++++ 3 files changed, 153 insertions(+), 40 deletions(-) create mode 100644 src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index a727dea4567..660b5fc556a 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -21,6 +21,8 @@ #include "common/DataTypes.hpp" #include "denseLinearAlgebra/common/layouts.hpp" +#include "LvArray/src/tensorOps.hpp" +#include "common/logger/Logger.hpp" #include @@ -30,7 +32,7 @@ namespace geos namespace denseLinearAlgebra { -namespace internal +namespace details { /** * @brief Solves a 2x2 linear system A * x = b. @@ -59,12 +61,14 @@ void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && LvArray::tensorOps::internal::checkSizes< 2 >( b ); LvArray::tensorOps::internal::checkSizes< 2 >( x ); - real64 const detA = A[0][0] * A[1][1] - A[0][1] * A[1][0]; + real64 const detA = LvArray::tensorOps::determinant< 2 >( A ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon; , "Singular system." ); + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon , "Singular system." ); - x[0] = (A[1][1] * b[0] - A[0][1] * b[1] ) / detA; - x[1] = (A[0][0] * b[1] - A[1][0] * b[0] ) / detA; + real64 const invA = 1.0 / detA; + + x[0] = ( A[1][1] * b[0] - A[0][1] * b[1] ) * invA; + x[1] = ( A[0][0] * b[1] - A[1][0] * b[0] ) * invA; } /** @@ -94,27 +98,27 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP LvArray::tensorOps::internal::checkSizes< 3 >( b ); LvArray::tensorOps::internal::checkSizes< 3 >( x ); - real64 const detA = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) + - A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]); + real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); + + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon, "Singular system." ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon; , "Singular system." ); + real64 const invA = 1.0 / detA; - real64 const detX0 = b[0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - - b[1] * (A[0][1] * A[2][2] - A[0][2] * A[2][1]) + - b[2] * (A[0][1] * A[1][2] - A[0][2] * A[1][1]); + real64 const detX0 = b[0] * ( A[1][1] * A[2][2] - A[2][1] * A[1][2] ) - + b[1] * ( A[0][1] * A[2][2] - A[0][2] * A[2][1] ) + + b[2] * ( A[0][1] * A[1][2] - A[0][2] * A[1][1] ); - real64 const detX1 = A[0][0] * (b[1] * A[2][2] - b[2] * A[2][1]) - - A[0][1] * (b[0] * A[2][2] - b[2] * A[2][0]) + - A[0][2] * (b[0] * A[1][2] - b[1] * A[1][0]); + real64 const detX1 = A[0][0] * ( b[1] * A[2][2] - b[2] * A[1][2] ) - + A[1][0] * ( b[0] * A[2][2] - b[2] * A[0][2] ) + + A[2][0] * ( b[0] * A[1][2] - b[1] * A[0][2] ); - real64 const detX2 = A[0][0] * (A[1][1] * b[2] - A[1][2] * b[1]) - - A[0][1] * (A[1][0] * b[2] - A[1][2] * b[0]) + - A[0][2] * (A[1][0] * b[1] - A[1][1] * b[0]); + real64 const detX2 = A[0][0] * ( A[1][1] * b[2] - A[2][1] * b[1] ) - + A[1][0] * ( A[0][1] * b[2] - A[2][1] * b[0] ) + + A[2][0] * ( A[0][1] * b[1] - A[1][1] * b[0] ); - x[0] = detX0 / detA; - x[1] = detX1 / detA; - x[2] = detX2 / detA; + x[0] = detX0 * invA; + x[1] = detX1 * invA; + x[2] = detX2 * invA; } /** @@ -128,14 +132,14 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP * @tparam RHS_TYPE The type of the right-hand side vector `b`. * @tparam SOL_TYPE The type of the solution vector `x`. * @param[in] A The upper triangular matrix representing the coefficients of the system. - * @param[in,out] b The right-hand side vector. It is used to compute the solution, but the values are not preserved. + * @param[in] b The right-hand side vector. It is used to compute the solution. * @param[out] x The solution vector. The result of solving the system `Ax = b` using back substitution. */ template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > -void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE & b, SOL_TYPE && x ) +void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { for( std::ptrdiff_t i = N - 1; i >= 0; --i ) { @@ -172,9 +176,9 @@ inline void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); - internal::checkSizes< N, N >( matrix ); - internal::checkSizes< N >( b ); - internal::checkSizes< N >( x ); + LvArray::tensorOps::internal::checkSizes< N, N >( A ); + LvArray::tensorOps::internal::checkSizes< N >( b ); + LvArray::tensorOps::internal::checkSizes< N >( x ); // Step 1: Transform into an upper triangular matrix @@ -213,7 +217,7 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) } // Step 2: Backward substitution - solveUpperTriangularSystem< N >( A, b, std::forward< N >( x ) ) + solveUpperTriangularSystem< N >( A, b, std::forward< N >( x ) ); } }; // internal namespace @@ -232,8 +236,8 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) * @tparam MODIFY_MATRIX Boolean flag indicating whether the input matrix `A` and vector `b` should be modified. * If `true`, the matrix `A` and vector `b` are modified in place. If `false`, copies of * `A` and `b` are made, and the original data is left unchanged. - * @param[in] A The constant matrix representing the coefficients of the system. - * @param[in] b The constant right-hand side vector. + * @param[in] A The matrix representing the coefficients of the system. + * @param[in] b The right-hand side vector. * @param[out] x The solution vector. The result of solving the system `Ax = b`. */ template< std::ptrdiff_t N, @@ -247,30 +251,30 @@ void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); static_assert( N < 10, "N cannot be larger than 9" ); - internal::checkSizes< N, N >( A ); - internal::checkSizes< N >( b ); - internal::checkSizes< N >( x ); + LvArray::tensorOps::internal::checkSizes< N, N >( A ); + LvArray::tensorOps::internal::checkSizes< N >( b ); + LvArray::tensorOps::internal::checkSizes< N >( x ); if constexpr ( N == 2 ) { - internal::solveTwoByTwoSystem( A, b, std::forward< SOL_TYPE >( x ) ); + details::solveTwoByTwoSystem( A, b, std::forward< SOL_TYPE >( x ) ); } else if constexpr ( N == 3 ) { - internal::solveThreeByThreeSystem( A, b, std::forward< SOL_TYPE >( x ) ); + details::solveThreeByThreeSystem( A, b, std::forward< SOL_TYPE >( x ) ); } else { if constexpr ( MODIFY_MATRIX ) { - internal::solveGaussianElimination< N >( A, b, std::forward< SOL_TYPE >( x ) ); + details::solveGaussianElimination< N >( A, b, std::forward< SOL_TYPE >( x ) ); } else { - real64[N][N] A_copy{}; - real64[N] b_copy{}; + real64 A_copy[N][N]{}; + real64 b_copy[N]{}; - for( std::ptrdiff_t i=0; i < N; ++j ) + for( std::ptrdiff_t i=0; i < N; ++i ) { b_copy[i] = b[i]; for( std::ptrdiff_t j=0; j < N; ++j ) @@ -278,7 +282,7 @@ void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) A_copy[i][j] = A[i][j]; } } - internal::solveGaussianElimination< N >( A_copy, b_copy, std::forward< SOL_TYPE >( x ) ); + details::solveGaussianElimination< N >( A_copy, b_copy, std::forward< SOL_TYPE >( x ) ); } } } diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt b/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt index 114dee90aac..5334abc1acb 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt +++ b/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt @@ -1,6 +1,7 @@ set( serial_tests testBlasLapack.cpp - testSolveLinearSystem.cpp ) + testSolveLinearSystem.cpp + testDenseLASolvers.cpp ) set( dependencyList gtest denseLinearAlgebra ) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp new file mode 100644 index 00000000000..1cb0a7bd080 --- /dev/null +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -0,0 +1,108 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "denseLinearAlgebra/denseLASolvers.hpp" + +// TPL includes +#include + + +using namespace geos; + + +constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon; + +template< std::ptrdiff_t N > +struct Reference +{}; + +template<> +struct Reference< 2 > +{ + + static constexpr std::ptrdiff_t size() { return 2; } + + static constexpr real64 A[2][2] = { { 2.0, -1.0}, + { 0.0, 1.0 } }; + + static constexpr real64 rhs[2] = {0.0, 1.0}; + + static constexpr real64 solution[2] = {0.5, 1.0}; +}; + +template<> +struct Reference< 3 > +{ + static constexpr std::ptrdiff_t size() { return 3; } + + static constexpr real64 A[3][3] = { { 4.0, 3.0, -2.0}, + { 2.0, -1.0, 5.0 }, + { -1.0, 3.0, -2.0 } }; + + static constexpr real64 rhs[3] = {-8.0, 19.0, -13.0}; + + static constexpr real64 solution[3] = {1.0, -2.0, 3.0}; +}; + +template< typename MATRIX_TYPE, typename VECTOR_TYPE, typename REFERENCE > +void assemble( MATRIX_TYPE & A, VECTOR_TYPE & rhs ) +{ + for( std::ptrdiff_t i=0; i < REFERENCE::size(); ++i ) + { + rhs[i] = REFERENCE::rhs[i]; + for( std::ptrdiff_t j=0; j < REFERENCE::size(); ++j ) + { + A[i][j] = REFERENCE::A[i][j]; + } + } +} + +template< typename REFERENCE > +void test_denseLASolve() +{ + constexpr std::ptrdiff_t N = REFERENCE::size(); + + real64 A[N][N]; + real64 b[N]; + real64 sol[N]; + + assemble< decltype(A), decltype(b), REFERENCE >( A, b ); + + denseLinearAlgebra::solve< N >( A, b, sol ); + + for( std::ptrdiff_t i = 0; i < N; ++i ) + { + EXPECT_NEAR( sol[i], + REFERENCE::solution[i], + machinePrecision ); + } +} + +TEST( denseLASolve, testTwoByTwo ) +{ + test_denseLASolve< Reference< 2 > >(); +} + +TEST( denseLASolve, testThreeByThree ) +{ + test_denseLASolve< Reference< 3 > >(); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} \ No newline at end of file From 05da4fd8414b69fb9325d4fb645858523a7b33c3 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 19 Aug 2024 21:57:22 -0700 Subject: [PATCH 07/21] now make system random. --- .../denseLinearAlgebra/denseLASolvers.hpp | 10 +-- .../unitTests/testDenseLASolvers.cpp | 84 ++++++++----------- 2 files changed, 41 insertions(+), 53 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 660b5fc556a..e62bfa5626a 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -61,9 +61,9 @@ void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && LvArray::tensorOps::internal::checkSizes< 2 >( b ); LvArray::tensorOps::internal::checkSizes< 2 >( x ); - real64 const detA = LvArray::tensorOps::determinant< 2 >( A ); + real64 const detA = LvArray::tensorOps::determinant< 2 >( A ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon , "Singular system." ); + GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon, "Singular system." ); real64 const invA = 1.0 / detA; @@ -98,7 +98,7 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP LvArray::tensorOps::internal::checkSizes< 3 >( b ); LvArray::tensorOps::internal::checkSizes< 3 >( x ); - real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); + real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon, "Singular system." ); @@ -110,7 +110,7 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP real64 const detX1 = A[0][0] * ( b[1] * A[2][2] - b[2] * A[1][2] ) - A[1][0] * ( b[0] * A[2][2] - b[2] * A[0][2] ) + - A[2][0] * ( b[0] * A[1][2] - b[1] * A[0][2] ); + A[2][0] * ( b[0] * A[1][2] - b[1] * A[0][2] ); real64 const detX2 = A[0][0] * ( A[1][1] * b[2] - A[2][1] * b[1] ) - A[1][0] * ( A[0][1] * b[2] - A[2][1] * b[0] ) + @@ -234,7 +234,7 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) * @tparam RHS_TYPE The type of the right-hand side vector `b`. * @tparam SOL_TYPE The type of the solution vector `x`. * @tparam MODIFY_MATRIX Boolean flag indicating whether the input matrix `A` and vector `b` should be modified. - * If `true`, the matrix `A` and vector `b` are modified in place. If `false`, copies of + * If `true`, the matrix `A` and vector `b` are modified in place. If `false`, copies of * `A` and `b` are made, and the original data is left unchanged. * @param[in] A The matrix representing the coefficients of the system. * @param[in] b The right-hand side vector. diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index 1cb0a7bd080..e4137e073d0 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -18,6 +18,8 @@ // TPL includes #include +#include + using namespace geos; @@ -25,79 +27,65 @@ using namespace geos; constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon; template< std::ptrdiff_t N > -struct Reference -{}; - -template<> -struct Reference< 2 > +struct LinearSystem { + static constexpr std::ptrdiff_t size() { return N; } - static constexpr std::ptrdiff_t size() { return 2; } - - static constexpr real64 A[2][2] = { { 2.0, -1.0}, - { 0.0, 1.0 } }; - - static constexpr real64 rhs[2] = {0.0, 1.0}; - - static constexpr real64 solution[2] = {0.5, 1.0}; -}; - -template<> -struct Reference< 3 > -{ - static constexpr std::ptrdiff_t size() { return 3; } - - static constexpr real64 A[3][3] = { { 4.0, 3.0, -2.0}, - { 2.0, -1.0, 5.0 }, - { -1.0, 3.0, -2.0 } }; - - static constexpr real64 rhs[3] = {-8.0, 19.0, -13.0}; - - static constexpr real64 solution[3] = {1.0, -2.0, 3.0}; -}; - -template< typename MATRIX_TYPE, typename VECTOR_TYPE, typename REFERENCE > -void assemble( MATRIX_TYPE & A, VECTOR_TYPE & rhs ) -{ - for( std::ptrdiff_t i=0; i < REFERENCE::size(); ++i ) + LinearSystem( short int seed ) { - rhs[i] = REFERENCE::rhs[i]; - for( std::ptrdiff_t j=0; j < REFERENCE::size(); ++j ) + std::mt19937 generator( seed ); + std::uniform_real_distribution< real64 > distribution( -10.0, 10.0 ); + for( ptrdiff_t i=0; i + real64 matrix[N][N]; + real64 solution[N]; + real64 rhs[N]; +}; + +template< typename LINEAR_SYSTEM > void test_denseLASolve() { - constexpr std::ptrdiff_t N = REFERENCE::size(); + constexpr std::ptrdiff_t N = LINEAR_SYSTEM::size(); - real64 A[N][N]; - real64 b[N]; + LINEAR_SYSTEM LS(2024); real64 sol[N]; - assemble< decltype(A), decltype(b), REFERENCE >( A, b ); - - denseLinearAlgebra::solve< N >( A, b, sol ); + denseLinearAlgebra::solve< N >( LS.matrix, LS.rhs, sol ); for( std::ptrdiff_t i = 0; i < N; ++i ) { EXPECT_NEAR( sol[i], - REFERENCE::solution[i], + LS.solution[i], machinePrecision ); } } TEST( denseLASolve, testTwoByTwo ) { - test_denseLASolve< Reference< 2 > >(); + test_denseLASolve< LinearSystem< 2 > >(); } TEST( denseLASolve, testThreeByThree ) { - test_denseLASolve< Reference< 3 > >(); + test_denseLASolve< LinearSystem< 3 > >(); } int main( int argc, char * * argv ) @@ -105,4 +93,4 @@ int main( int argc, char * * argv ) ::testing::InitGoogleTest( &argc, argv ); int const result = RUN_ALL_TESTS(); return result; -} \ No newline at end of file +} From b1d9fe05fc0131db7644b2263592379ff43b27e4 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 10:28:37 -0700 Subject: [PATCH 08/21] unitTests work fine now. --- .../denseLinearAlgebra/denseLASolvers.hpp | 43 ++++-- .../unitTests/testDenseLASolvers.cpp | 139 +++++++++++++++--- 2 files changed, 145 insertions(+), 37 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index e62bfa5626a..31bb7697803 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -34,6 +34,9 @@ namespace denseLinearAlgebra namespace details { + +constexpr real64 singularMatrixTolerance = 1e2*LvArray::NumericLimits< real64 >::epsilon; + /** * @brief Solves a 2x2 linear system A * x = b. * @@ -49,13 +52,14 @@ namespace details * @param[in] A The 2x2 matrix representing the system of equations. Must have size 2x2. * @param[in] b The 2-element vector representing the right-hand side of the equation. * @param[out] x The 2-element vector that will store the solution to the system. + * @return bool whether the solve succeeded or not. */ template< typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +bool solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { LvArray::tensorOps::internal::checkSizes< 2, 2 >( A ); LvArray::tensorOps::internal::checkSizes< 2 >( b ); @@ -63,12 +67,15 @@ void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && real64 const detA = LvArray::tensorOps::determinant< 2 >( A ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon, "Singular system." ); + if( LvArray::math::abs( detA ) < singularMatrixTolerance ) + return false; real64 const invA = 1.0 / detA; x[0] = ( A[1][1] * b[0] - A[0][1] * b[1] ) * invA; x[1] = ( A[0][0] * b[1] - A[1][0] * b[0] ) * invA; + + return true; } /** @@ -86,13 +93,14 @@ void solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && * @param[in] A The 3x3 matrix representing the system of equations. Must have size 3x3. * @param[in] b The 3-element vector representing the right-hand side of the equation. * @param[out] x The 3-element vector that will store the solution to the system. + * @return bool whether the solve succeeded or not. */ template< typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +bool solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { LvArray::tensorOps::internal::checkSizes< 3, 3 >( A ); LvArray::tensorOps::internal::checkSizes< 3 >( b ); @@ -100,7 +108,10 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( detA ), LvArray::NumericLimits< real64 >::epsilon, "Singular system." ); + std::cout << "detA: " << detA << std::endl; + + if( LvArray::math::abs( detA ) < singularMatrixTolerance ) + return false; real64 const invA = 1.0 / detA; @@ -119,6 +130,8 @@ void solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP x[0] = detX0 * invA; x[1] = detX1 * invA; x[2] = detX2 * invA; + + return true; } /** @@ -166,6 +179,7 @@ void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_ * @param[in,out] A The matrix to be transformed into an upper triangular matrix. Modified in place. * @param[in,out] b The right-hand side vector. Modified in place to reflect the transformed system. * @param[out] x The solution vector. The result of solving the system `Ax = b`. + * @return bool whether the solve succeeded or not. */ template< std::ptrdiff_t N, typename MATRIX_TYPE, @@ -173,7 +187,7 @@ template< std::ptrdiff_t N, typename SOL_TYPE > GEOS_HOST_DEVICE inline -void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); LvArray::tensorOps::internal::checkSizes< N, N >( A ); @@ -202,7 +216,8 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) } std::swap( b[i], b[max_row] ); - GEOS_ERROR_IF_LT_MSG( LvArray::math::abs( A[i][i] ), LvArray::NumericLimits< real64 >::epsilon, "Singular matrix." ); + if( LvArray::math::abs( A[i][i] ) < singularMatrixTolerance ) + return false; // 1.c Eliminate entries below the pivot for( std::ptrdiff_t k = i + 1; k < N; ++k ) @@ -217,7 +232,9 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) } // Step 2: Backward substitution - solveUpperTriangularSystem< N >( A, b, std::forward< N >( x ) ); + solveUpperTriangularSystem< N >( A, b, std::forward< SOL_TYPE >( x ) ); + + return true; } }; // internal namespace @@ -239,6 +256,7 @@ void solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) * @param[in] A The matrix representing the coefficients of the system. * @param[in] b The right-hand side vector. * @param[out] x The solution vector. The result of solving the system `Ax = b`. + * @return bool whether the solve succeeded or not. */ template< std::ptrdiff_t N, typename MATRIX_TYPE, @@ -247,7 +265,7 @@ template< std::ptrdiff_t N, bool MODIFY_MATRIX = true > GEOS_HOST_DEVICE inline -void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +bool solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); static_assert( N < 10, "N cannot be larger than 9" ); @@ -257,17 +275,17 @@ void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) if constexpr ( N == 2 ) { - details::solveTwoByTwoSystem( A, b, std::forward< SOL_TYPE >( x ) ); + return details::solveTwoByTwoSystem( A, b, std::forward< SOL_TYPE >( x ) ); } else if constexpr ( N == 3 ) { - details::solveThreeByThreeSystem( A, b, std::forward< SOL_TYPE >( x ) ); + return details::solveThreeByThreeSystem( A, b, std::forward< SOL_TYPE >( x ) ); } else { if constexpr ( MODIFY_MATRIX ) { - details::solveGaussianElimination< N >( A, b, std::forward< SOL_TYPE >( x ) ); + return details::solveGaussianElimination< N >( A, b, std::forward< SOL_TYPE >( x ) ); } else { @@ -282,11 +300,10 @@ void solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) A_copy[i][j] = A[i][j]; } } - details::solveGaussianElimination< N >( A_copy, b_copy, std::forward< SOL_TYPE >( x ) ); + return details::solveGaussianElimination< N >( A_copy, b_copy, std::forward< SOL_TYPE >( x ) ); } } } - }; }; diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index e4137e073d0..59763c33d0c 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -20,27 +20,43 @@ #include - using namespace geos; constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon; + template< std::ptrdiff_t N > struct LinearSystem { - static constexpr std::ptrdiff_t size() { return N; } +public: + real64 matrix[N][N]; + real64 solution[N]; + real64 rhs[N]; +}; - LinearSystem( short int seed ) +template< std::ptrdiff_t N > +struct InvertibleLinearSystem : public LinearSystem< N > +{ + using LinearSystem< N >::matrix; + using LinearSystem< N >::solution; + using LinearSystem< N >::rhs; + InvertibleLinearSystem( short int seed ) { - std::mt19937 generator( seed ); + std::mt19937 generator( seed ); std::uniform_real_distribution< real64 > distribution( -10.0, 10.0 ); + std::uniform_real_distribution< real64 > perturbation( -20.0, 20.0 ); for( ptrdiff_t i=0; i +struct SingularLinearSystem : public LinearSystem< N > +{ + using LinearSystem< N >::matrix; + using LinearSystem< N >::solution; + using LinearSystem< N >::rhs; + SingularLinearSystem( short int seed ) + { + std::mt19937 generator( seed ); + std::uniform_real_distribution< real64 > distribution( -10.0, 10.0 ); + for( ptrdiff_t i=0; i -void test_denseLASolve() +template< typename N > +class DenseLinearSolverTest : public ::testing::Test { - constexpr std::ptrdiff_t N = LINEAR_SYSTEM::size(); +public: - LINEAR_SYSTEM LS(2024); - real64 sol[N]; + static constexpr std::ptrdiff_t size = N::value; - denseLinearAlgebra::solve< N >( LS.matrix, LS.rhs, sol ); + DenseLinearSolverTest() = default; + ~DenseLinearSolverTest() override = default; - for( std::ptrdiff_t i = 0; i < N; ++i ) + void test_solve() { - EXPECT_NEAR( sol[i], - LS.solution[i], - machinePrecision ); + InvertibleLinearSystem< size > LS( 2024 ); + real64 sol[size]; + + bool const success = denseLinearAlgebra::solve< size >( LS.matrix, LS.rhs, sol ); + + EXPECT_TRUE( success ); + + for( std::ptrdiff_t i = 0; i < size; ++i ) + { + EXPECT_NEAR( sol[i], + LS.solution[i], + machinePrecision ); + } } -} + void test_singularSystem() + { + SingularLinearSystem< size > LS( 2024 ); + real64 sol[size]; + + bool const success = denseLinearAlgebra::solve< size >( LS.matrix, LS.rhs, sol ); + + EXPECT_FALSE( success ); + } + +}; -TEST( denseLASolve, testTwoByTwo ) +using Dimensions = ::testing::Types< std::integral_constant< std::ptrdiff_t, 2 >, + std::integral_constant< std::ptrdiff_t, 3 >, + std::integral_constant< std::ptrdiff_t, 4 >, + std::integral_constant< std::ptrdiff_t, 5 >, + std::integral_constant< std::ptrdiff_t, 6 >, + std::integral_constant< std::ptrdiff_t, 7 >, + std::integral_constant< std::ptrdiff_t, 8 >, + std::integral_constant< std::ptrdiff_t, 9 > >; + +class NameGenerator { - test_denseLASolve< LinearSystem< 2 > >(); -} +public: + template< typename T > + static std::string GetName( int ) + { + if constexpr (T::value == 2) return "TwoByTwo"; + if constexpr (T::value == 3) return "ThreeByThree"; + if constexpr (T::value == 4) return "FourByFour"; + if constexpr (T::value == 5) return "FiveByFive"; + if constexpr (T::value == 6) return "SixBySix"; + if constexpr (T::value == 7) return "SevenBySeven"; + if constexpr (T::value == 8) return "EightByEight"; + if constexpr (T::value == 9) return "NineByNine"; + } +}; + +TYPED_TEST_SUITE( DenseLinearSolverTest, Dimensions, NameGenerator ); -TEST( denseLASolve, testThreeByThree ) +TYPED_TEST( DenseLinearSolverTest, testDenseLA ) { - test_denseLASolve< LinearSystem< 3 > >(); + this->test_solve(); + this->test_singularSystem(); } int main( int argc, char * * argv ) From fcf8d6fe6bd35078cc47347d6bdc6c5020b46118 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 10:30:42 -0700 Subject: [PATCH 09/21] clean up --- src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 31bb7697803..c0f0e01d05c 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -108,8 +108,6 @@ bool solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); - std::cout << "detA: " << detA << std::endl; - if( LvArray::math::abs( detA ) < singularMatrixTolerance ) return false; @@ -194,7 +192,6 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) LvArray::tensorOps::internal::checkSizes< N >( b ); LvArray::tensorOps::internal::checkSizes< N >( x ); - // Step 1: Transform into an upper triangular matrix // 1.a. Find the pivot From ba99762812b6412f13e4c7e166046506cb000ef3 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 10:44:51 -0700 Subject: [PATCH 10/21] small fixes. --- .../denseLinearAlgebra/unitTests/testDenseLASolvers.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index 59763c33d0c..cd4889049ce 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -52,12 +52,9 @@ struct InvertibleLinearSystem : public LinearSystem< N > for( ptrdiff_t j=0; j Date: Tue, 20 Aug 2024 10:48:34 -0700 Subject: [PATCH 11/21] Update src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp --- src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index c0f0e01d05c..296b3d4b7d7 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -234,7 +234,7 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) return true; } -}; // internal namespace +}; // details namespace /** * @brief Solves a linear system using the most appropriate method based on the size of the system. From b48de388e22edc2eebba6829450658519d89d75a Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 15:39:45 -0700 Subject: [PATCH 12/21] add test utils. --- .../unitTests/testDenseLASolvers.cpp | 59 ++++++++++++++----- .../unitTests/testUtils.hpp | 34 +++++++++++ 2 files changed, 79 insertions(+), 14 deletions(-) create mode 100644 src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index cd4889049ce..c08e1ea3445 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -14,17 +14,24 @@ */ #include "denseLinearAlgebra/denseLASolvers.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "testUtils.hpp" // TPL includes #include #include -using namespace geos; +namespace geos +{ +namespace denseLinearAlgebra +{ -constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon; +namespace testing +{ +constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon; template< std::ptrdiff_t N > struct LinearSystem @@ -117,27 +124,44 @@ class DenseLinearSolverTest : public ::testing::Test void test_solve() { InvertibleLinearSystem< size > LS( 2024 ); - real64 sol[size]; - bool const success = denseLinearAlgebra::solve< size >( LS.matrix, LS.rhs, sol ); + forAll< parallelDevicePolicy<> >( 1, GEOS_HOST_DEVICE [=]( int ) + { + real64 sol[size]; + real64 matrix[size][size]; + real64 rhs[size]; - EXPECT_TRUE( success ); + LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); + LvArray::tensorOps::copy< size >( rhs, LS.rhs ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - for( std::ptrdiff_t i = 0; i < size; ++i ) - { - EXPECT_NEAR( sol[i], - LS.solution[i], - machinePrecision ); - } + EXPECT_TRUE( success ); + + for( std::ptrdiff_t i = 0; i < size; ++i ) + { + PORTABLE_EXPECT_NEAR( sol[i], + LS.solution[i], + machinePrecision ); + } + } ); } void test_singularSystem() { SingularLinearSystem< size > LS( 2024 ); - real64 sol[size]; - bool const success = denseLinearAlgebra::solve< size >( LS.matrix, LS.rhs, sol ); + forAll< parallelDevicePolicy<> >( 1, GEOS_HOST_DEVICE [=]( int ) + { + + real64 sol[size]; + real64 matrix[size][size]; + real64 rhs[size]; + + LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); + LvArray::tensorOps::copy< size >( rhs, LS.rhs ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - EXPECT_FALSE( success ); + EXPECT_FALSE( success ); + } ); } }; @@ -182,3 +206,10 @@ int main( int argc, char * * argv ) int const result = RUN_ALL_TESTS(); return result; } + + +} // testing + +} // denseLinearAlgebra + +} // namespace geos \ No newline at end of file diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp new file mode 100644 index 00000000000..3da0632f7c5 --- /dev/null +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp @@ -0,0 +1,34 @@ + +#include "common/GeosxMacros.hpp" + +// TPL includes +#include + + +namespace geos +{ +namespace denseLinearAlgebra +{ +namespace testing +{ + + +#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__) +#define PORTABLE_EXPECT_EQ( L, R ) GEOS_ERROR_IF_NE( L, R ) +#define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) LVARRAY_ERROR_IF_GE_MSG( LvArray::math::abs( ( L ) -( R ) ), EPSILON, \ + STRINGIZE( L ) " = " << ( L ) << "\n" << STRINGIZE( R ) " = " << ( R ) ); +#define PORTABLE_EXPECT_TRUE( value ) GEOS_ERROR_IF( value, "should be true" ); +#define PORTABLE_EXPECT_FALSE( value ) GEOS_ERROR_IF( !value, "should be false" ); +#else +#define PORTABLE_EXPECT_EQ( L, R ) EXPECT_EQ( L, R ) +#define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) EXPECT_LE( LvArray::math::abs( ( L ) -( R ) ), EPSILON ) << \ + STRINGIZE( L ) " = " << ( L ) << "\n" << STRINGIZE( R ) " = " << ( R ); +#define PORTABLE_EXPECT_TRUE( value ) EXPECT_TRUE( value ) +#define PORTABLE_EXPECT_FALSE( value ) EXPECT_FALSE( value ) +#endif + +} //namespace testing + +} // namespace denseLinearAlgebra + +} // namespace geos From 987caf7611166ed8003ab4f53501ffeeb44bbfb1 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 16:35:58 -0700 Subject: [PATCH 13/21] use GEOS_DEVICE_COMPILE macro instead --- .../unitTests/testDenseLASolvers.cpp | 60 +++++++++---------- .../unitTests/testUtils.hpp | 8 +-- 2 files changed, 32 insertions(+), 36 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index c08e1ea3445..2c3dc24d867 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -24,10 +24,8 @@ namespace geos { - namespace denseLinearAlgebra { - namespace testing { @@ -75,7 +73,6 @@ struct InvertibleLinearSystem : public LinearSystem< N > } } }; - template< std::ptrdiff_t N > struct SingularLinearSystem : public LinearSystem< N > { @@ -126,42 +123,42 @@ class DenseLinearSolverTest : public ::testing::Test InvertibleLinearSystem< size > LS( 2024 ); forAll< parallelDevicePolicy<> >( 1, GEOS_HOST_DEVICE [=]( int ) - { - real64 sol[size]; - real64 matrix[size][size]; - real64 rhs[size]; - - LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); - LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - - EXPECT_TRUE( success ); - - for( std::ptrdiff_t i = 0; i < size; ++i ) - { - PORTABLE_EXPECT_NEAR( sol[i], - LS.solution[i], - machinePrecision ); - } - } ); + { + real64 sol[size]; + real64 matrix[size][size]; + real64 rhs[size]; + + LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); + LvArray::tensorOps::copy< size >( rhs, LS.rhs ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + + EXPECT_TRUE( success ); + + for( std::ptrdiff_t i = 0; i < size; ++i ) + { + PORTABLE_EXPECT_NEAR( sol[i], + LS.solution[i], + machinePrecision ); + } + } ); } void test_singularSystem() { SingularLinearSystem< size > LS( 2024 ); forAll< parallelDevicePolicy<> >( 1, GEOS_HOST_DEVICE [=]( int ) - { + { - real64 sol[size]; - real64 matrix[size][size]; - real64 rhs[size]; + real64 sol[size]; + real64 matrix[size][size]; + real64 rhs[size]; - LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); - LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); + LvArray::tensorOps::copy< size >( rhs, LS.rhs ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - EXPECT_FALSE( success ); - } ); + EXPECT_FALSE( success ); + } ); } }; @@ -207,9 +204,8 @@ int main( int argc, char * * argv ) return result; } - } // testing } // denseLinearAlgebra -} // namespace geos \ No newline at end of file +} // namespace geos diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp index 3da0632f7c5..e28f68ac555 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp @@ -8,23 +8,23 @@ namespace geos { namespace denseLinearAlgebra -{ +{ namespace testing { -#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__) +#if defined(GEOS_DEVICE_COMPILE) #define PORTABLE_EXPECT_EQ( L, R ) GEOS_ERROR_IF_NE( L, R ) #define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) LVARRAY_ERROR_IF_GE_MSG( LvArray::math::abs( ( L ) -( R ) ), EPSILON, \ STRINGIZE( L ) " = " << ( L ) << "\n" << STRINGIZE( R ) " = " << ( R ) ); #define PORTABLE_EXPECT_TRUE( value ) GEOS_ERROR_IF( value, "should be true" ); -#define PORTABLE_EXPECT_FALSE( value ) GEOS_ERROR_IF( !value, "should be false" ); +#define PORTABLE_EXPECT_FALSE( value ) GEOS_ERROR_IF( !value, "should be false" ); #else #define PORTABLE_EXPECT_EQ( L, R ) EXPECT_EQ( L, R ) #define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) EXPECT_LE( LvArray::math::abs( ( L ) -( R ) ), EPSILON ) << \ STRINGIZE( L ) " = " << ( L ) << "\n" << STRINGIZE( R ) " = " << ( R ); #define PORTABLE_EXPECT_TRUE( value ) EXPECT_TRUE( value ) -#define PORTABLE_EXPECT_FALSE( value ) EXPECT_FALSE( value ) +#define PORTABLE_EXPECT_FALSE( value ) EXPECT_FALSE( value ) #endif } //namespace testing From 1d2e95c98c5f4e2ffd096bc8568e969f0181200f Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 16:38:30 -0700 Subject: [PATCH 14/21] remove commas. --- src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 296b3d4b7d7..8b4143f439a 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -234,7 +234,7 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) return true; } -}; // details namespace +} // details namespace /** * @brief Solves a linear system using the most appropriate method based on the size of the system. @@ -301,9 +301,10 @@ bool solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) } } } -}; -}; +} // denseLinearAlgebra + +} // geos #endif /*GEOS_DENSELINEARALGEBRA_DENSELASOLVERS_HPP_*/ From 691afbbf1617bfcbed63f8d7f27ea260fdfbca4d Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 16:40:57 -0700 Subject: [PATCH 15/21] modify tolerance. --- .../denseLinearAlgebra/unitTests/testDenseLASolvers.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index 2c3dc24d867..bd4f8df752f 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -29,7 +29,7 @@ namespace denseLinearAlgebra namespace testing { -constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon; +constexpr real64 machinePrecision = 1.0e3 * LvArray::NumericLimits< real64 >::epsilon; template< std::ptrdiff_t N > struct LinearSystem From 9c4594aedf647955c096a81965ad7bafed492282 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 16:55:54 -0700 Subject: [PATCH 16/21] uncrustify --- .../unitTests/testDenseLASolvers.cpp | 58 +++++++++---------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index bd4f8df752f..d98aa85c9ae 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -122,43 +122,43 @@ class DenseLinearSolverTest : public ::testing::Test { InvertibleLinearSystem< size > LS( 2024 ); - forAll< parallelDevicePolicy<> >( 1, GEOS_HOST_DEVICE [=]( int ) - { - real64 sol[size]; - real64 matrix[size][size]; - real64 rhs[size]; - - LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); - LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - - EXPECT_TRUE( success ); - - for( std::ptrdiff_t i = 0; i < size; ++i ) - { - PORTABLE_EXPECT_NEAR( sol[i], - LS.solution[i], - machinePrecision ); - } - } ); + forAll< parallelDevicePolicy<> >( 1, [=] GEOS_HOST_DEVICE ( int ) + { + real64 sol[size]; + real64 matrix[size][size]; + real64 rhs[size]; + + LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); + LvArray::tensorOps::copy< size >( rhs, LS.rhs ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + + EXPECT_TRUE( success ); + + for( std::ptrdiff_t i = 0; i < size; ++i ) + { + PORTABLE_EXPECT_NEAR( sol[i], + LS.solution[i], + machinePrecision ); + } + } ); } void test_singularSystem() { SingularLinearSystem< size > LS( 2024 ); - forAll< parallelDevicePolicy<> >( 1, GEOS_HOST_DEVICE [=]( int ) - { + forAll< parallelDevicePolicy<> >( 1, [=] GEOS_HOST_DEVICE ( int ) + { - real64 sol[size]; - real64 matrix[size][size]; - real64 rhs[size]; + real64 sol[size]; + real64 matrix[size][size]; + real64 rhs[size]; - LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); - LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); + LvArray::tensorOps::copy< size >( rhs, LS.rhs ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - EXPECT_FALSE( success ); - } ); + EXPECT_FALSE( success ); + } ); } }; From a710414242ce39697c4ba577a0bcd21b4467a550 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 17:05:06 -0700 Subject: [PATCH 17/21] portable bool check. --- .../denseLinearAlgebra/unitTests/testDenseLASolvers.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index d98aa85c9ae..83e7684e249 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -132,7 +132,7 @@ class DenseLinearSolverTest : public ::testing::Test LvArray::tensorOps::copy< size >( rhs, LS.rhs ); bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - EXPECT_TRUE( success ); + PORTABLE_EXPECT_TRUE( success ); for( std::ptrdiff_t i = 0; i < size; ++i ) { @@ -157,7 +157,7 @@ class DenseLinearSolverTest : public ::testing::Test LvArray::tensorOps::copy< size >( rhs, LS.rhs ); bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); - EXPECT_FALSE( success ); + PORTABLE_EXPECT_FALSE( success ); } ); } From 66b1bde23db2e1f140d117c94b5f5c9af097c02e Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 17:09:46 -0700 Subject: [PATCH 18/21] forgotten GEOS_HOST_DEVICE --- src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 8b4143f439a..adc3cd2e4d1 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -150,6 +150,8 @@ template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > +GEOS_HOST_DEVICE +inline void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { for( std::ptrdiff_t i = N - 1; i >= 0; --i ) From 0f9d1284dab7440b2404a3c9885422ceaff0582a Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 20 Aug 2024 17:34:38 -0700 Subject: [PATCH 19/21] remove std::swap --- .../denseLinearAlgebra/denseLASolvers.hpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index adc3cd2e4d1..3594a374dc4 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -202,7 +202,7 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) std::ptrdiff_t max_row = i; for( std::ptrdiff_t k = i + 1; k < N; ++k ) { - if( std::abs( A[k][i] ) > std::abs( A[max_row][i] )) + if( LvArray::math::abs( A[k][i] ) > LvArray::math::abs( A[max_row][i] )) { max_row = k; } @@ -211,9 +211,16 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) // 1.b. Swap rows for( std::ptrdiff_t k = i; k < N; ++k ) { - std::swap( A[i][k], A[max_row][k] ); + // std::swap( A[i][k], A[max_row][k] ); + real64 const temp = A[max_row][k]; + A[max_row][k] = A[i][k]; + A[i][k] = temp; } - std::swap( b[i], b[max_row] ); + // std::swap( b[i], b[max_row] ); cannot be done on device + real64 const temp = b[i]; + b[i] = b[max_row]; + b[max_row] = temp; + if( LvArray::math::abs( A[i][i] ) < singularMatrixTolerance ) return false; From 71d48a383d5605913cbc2ff2034e02ece525c84f Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Wed, 21 Aug 2024 10:54:54 -0700 Subject: [PATCH 20/21] use int instead of bool to avoid cuda error. --- .../denseLinearAlgebra/denseLASolvers.hpp | 36 ++++++------ .../unitTests/testDenseLASolvers.cpp | 56 ++++++++++++++++--- .../unitTests/testUtils.hpp | 4 +- 3 files changed, 67 insertions(+), 29 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index 3594a374dc4..b0bf28fb7ca 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -52,14 +52,14 @@ constexpr real64 singularMatrixTolerance = 1e2*LvArray::NumericLimits< real64 >: * @param[in] A The 2x2 matrix representing the system of equations. Must have size 2x2. * @param[in] b The 2-element vector representing the right-hand side of the equation. * @param[out] x The 2-element vector that will store the solution to the system. - * @return bool whether the solve succeeded or not. + * @return int that sepcifies whether the solve succeeded (1) or not (0). */ template< typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -bool solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +int solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { LvArray::tensorOps::internal::checkSizes< 2, 2 >( A ); LvArray::tensorOps::internal::checkSizes< 2 >( b ); @@ -68,14 +68,14 @@ bool solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && real64 const detA = LvArray::tensorOps::determinant< 2 >( A ); if( LvArray::math::abs( detA ) < singularMatrixTolerance ) - return false; + return 0; real64 const invA = 1.0 / detA; x[0] = ( A[1][1] * b[0] - A[0][1] * b[1] ) * invA; x[1] = ( A[0][0] * b[1] - A[1][0] * b[0] ) * invA; - return true; + return 1; } /** @@ -93,14 +93,14 @@ bool solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && * @param[in] A The 3x3 matrix representing the system of equations. Must have size 3x3. * @param[in] b The 3-element vector representing the right-hand side of the equation. * @param[out] x The 3-element vector that will store the solution to the system. - * @return bool whether the solve succeeded or not. + * @return int that sepcifies whether the solve succeeded (1) or not (0). */ template< typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -bool solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +int solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { LvArray::tensorOps::internal::checkSizes< 3, 3 >( A ); LvArray::tensorOps::internal::checkSizes< 3 >( b ); @@ -109,7 +109,7 @@ bool solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); if( LvArray::math::abs( detA ) < singularMatrixTolerance ) - return false; + return 0; real64 const invA = 1.0 / detA; @@ -129,7 +129,7 @@ bool solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYP x[1] = detX1 * invA; x[2] = detX2 * invA; - return true; + return 1; } /** @@ -179,7 +179,7 @@ void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_ * @param[in,out] A The matrix to be transformed into an upper triangular matrix. Modified in place. * @param[in,out] b The right-hand side vector. Modified in place to reflect the transformed system. * @param[out] x The solution vector. The result of solving the system `Ax = b`. - * @return bool whether the solve succeeded or not. + * @return int that sepcifies whether the solve succeeded (1) or not (0). */ template< std::ptrdiff_t N, typename MATRIX_TYPE, @@ -187,7 +187,7 @@ template< std::ptrdiff_t N, typename SOL_TYPE > GEOS_HOST_DEVICE inline -bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +int solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); LvArray::tensorOps::internal::checkSizes< N, N >( A ); @@ -220,10 +220,10 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) real64 const temp = b[i]; b[i] = b[max_row]; b[max_row] = temp; - + if( LvArray::math::abs( A[i][i] ) < singularMatrixTolerance ) - return false; + return 0; // 1.c Eliminate entries below the pivot for( std::ptrdiff_t k = i + 1; k < N; ++k ) @@ -240,7 +240,7 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) // Step 2: Backward substitution solveUpperTriangularSystem< N >( A, b, std::forward< SOL_TYPE >( x ) ); - return true; + return 1; } } // details namespace @@ -256,22 +256,22 @@ bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. * @tparam SOL_TYPE The type of the solution vector `x`. - * @tparam MODIFY_MATRIX Boolean flag indicating whether the input matrix `A` and vector `b` should be modified. - * If `true`, the matrix `A` and vector `b` are modified in place. If `false`, copies of + * @tparam MODIFY_MATRIX intean flag indicating whether the input matrix `A` and vector `b` should be modified. + * If `1`, the matrix `A` and vector `b` are modified in place. If `0`, copies of * `A` and `b` are made, and the original data is left unchanged. * @param[in] A The matrix representing the coefficients of the system. * @param[in] b The right-hand side vector. * @param[out] x The solution vector. The result of solving the system `Ax = b`. - * @return bool whether the solve succeeded or not. + * @return int that sepcifies whether the solve succeeded (1) or not (0). */ template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE, - bool MODIFY_MATRIX = true > + int MODIFY_MATRIX = 1 > GEOS_HOST_DEVICE inline -bool solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +int solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); static_assert( N < 10, "N cannot be larger than 9" ); diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index 83e7684e249..b0cd4c1eb4e 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -32,7 +32,7 @@ namespace testing constexpr real64 machinePrecision = 1.0e3 * LvArray::NumericLimits< real64 >::epsilon; template< std::ptrdiff_t N > -struct LinearSystem +class LinearSystem { public: real64 matrix[N][N]; @@ -41,11 +41,28 @@ struct LinearSystem }; template< std::ptrdiff_t N > -struct InvertibleLinearSystem : public LinearSystem< N > +class InvertibleLinearSystem : public LinearSystem< N > { +public: using LinearSystem< N >::matrix; using LinearSystem< N >::solution; using LinearSystem< N >::rhs; + + /// Default copy constructor + InvertibleLinearSystem( InvertibleLinearSystem const & ) = default; + + /// Default move constructor + InvertibleLinearSystem( InvertibleLinearSystem && ) = default; + + /// Deleted default constructor + InvertibleLinearSystem() = delete; + + /// Deleted copy assignment operator + InvertibleLinearSystem & operator=( InvertibleLinearSystem const & ) = delete; + + /// Deleted move assignment operator + InvertibleLinearSystem & operator=( InvertibleLinearSystem && ) = delete; + InvertibleLinearSystem( short int seed ) { std::mt19937 generator( seed ); @@ -73,12 +90,30 @@ struct InvertibleLinearSystem : public LinearSystem< N > } } }; + template< std::ptrdiff_t N > -struct SingularLinearSystem : public LinearSystem< N > +class SingularLinearSystem : public LinearSystem< N > { +public: using LinearSystem< N >::matrix; using LinearSystem< N >::solution; using LinearSystem< N >::rhs; + + /// Default copy constructor + SingularLinearSystem( SingularLinearSystem const & ) = default; + + /// Default move constructor + SingularLinearSystem( SingularLinearSystem && ) = default; + + /// Deleted default constructor + SingularLinearSystem() = delete; + + /// Deleted copy assignment operator + SingularLinearSystem & operator=( SingularLinearSystem const & ) = delete; + + /// Deleted move assignment operator + SingularLinearSystem & operator=( SingularLinearSystem && ) = delete; + SingularLinearSystem( short int seed ) { std::mt19937 generator( seed ); @@ -122,15 +157,17 @@ class DenseLinearSolverTest : public ::testing::Test { InvertibleLinearSystem< size > LS( 2024 ); + // GEOS_UNUSED_VAR(LS); + forAll< parallelDevicePolicy<> >( 1, [=] GEOS_HOST_DEVICE ( int ) { - real64 sol[size]; - real64 matrix[size][size]; - real64 rhs[size]; + real64 sol[size]{}; + real64 matrix[size][size]{}; + real64 rhs[size]{}; LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + int const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); PORTABLE_EXPECT_TRUE( success ); @@ -155,7 +192,7 @@ class DenseLinearSolverTest : public ::testing::Test LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + int const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); PORTABLE_EXPECT_FALSE( success ); } ); @@ -172,6 +209,7 @@ using Dimensions = ::testing::Types< std::integral_constant< std::ptrdiff_t, 2 > std::integral_constant< std::ptrdiff_t, 8 >, std::integral_constant< std::ptrdiff_t, 9 > >; + class NameGenerator { public: @@ -194,7 +232,7 @@ TYPED_TEST_SUITE( DenseLinearSolverTest, Dimensions, NameGenerator ); TYPED_TEST( DenseLinearSolverTest, testDenseLA ) { this->test_solve(); - this->test_singularSystem(); + // this->test_singularSystem(); } int main( int argc, char * * argv ) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp index e28f68ac555..82a94ca86f3 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp @@ -17,8 +17,8 @@ namespace testing #define PORTABLE_EXPECT_EQ( L, R ) GEOS_ERROR_IF_NE( L, R ) #define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) LVARRAY_ERROR_IF_GE_MSG( LvArray::math::abs( ( L ) -( R ) ), EPSILON, \ STRINGIZE( L ) " = " << ( L ) << "\n" << STRINGIZE( R ) " = " << ( R ) ); -#define PORTABLE_EXPECT_TRUE( value ) GEOS_ERROR_IF( value, "should be true" ); -#define PORTABLE_EXPECT_FALSE( value ) GEOS_ERROR_IF( !value, "should be false" ); +#define PORTABLE_EXPECT_TRUE( value ) LVARRAY_ERROR_IF( value==0, "should be true" ) +#define PORTABLE_EXPECT_FALSE( value ) LVARRAY_ERROR_IF( value==1, "should be false" ) #else #define PORTABLE_EXPECT_EQ( L, R ) EXPECT_EQ( L, R ) #define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) EXPECT_LE( LvArray::math::abs( ( L ) -( R ) ), EPSILON ) << \ From b529c54261a6c79f11d53ea7dc17bc8978221416 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Wed, 21 Aug 2024 11:38:38 -0700 Subject: [PATCH 21/21] revert to bool. --- .../denseLinearAlgebra/denseLASolvers.hpp | 38 ++++++++--------- .../unitTests/testDenseLASolvers.cpp | 42 +++---------------- .../unitTests/testUtils.hpp | 4 +- 3 files changed, 27 insertions(+), 57 deletions(-) diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp index b0bf28fb7ca..d45ae76accb 100644 --- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp +++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp @@ -52,14 +52,14 @@ constexpr real64 singularMatrixTolerance = 1e2*LvArray::NumericLimits< real64 >: * @param[in] A The 2x2 matrix representing the system of equations. Must have size 2x2. * @param[in] b The 2-element vector representing the right-hand side of the equation. * @param[out] x The 2-element vector that will store the solution to the system. - * @return int that sepcifies whether the solve succeeded (1) or not (0). + * @return bool that sepcifies whether the solve succeeded (1) or not (0). */ template< typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -int solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +bool solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { LvArray::tensorOps::internal::checkSizes< 2, 2 >( A ); LvArray::tensorOps::internal::checkSizes< 2 >( b ); @@ -68,14 +68,14 @@ int solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && real64 const detA = LvArray::tensorOps::determinant< 2 >( A ); if( LvArray::math::abs( detA ) < singularMatrixTolerance ) - return 0; + return false; real64 const invA = 1.0 / detA; x[0] = ( A[1][1] * b[0] - A[0][1] * b[1] ) * invA; x[1] = ( A[0][0] * b[1] - A[1][0] * b[0] ) * invA; - return 1; + return true; } /** @@ -93,14 +93,14 @@ int solveTwoByTwoSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && * @param[in] A The 3x3 matrix representing the system of equations. Must have size 3x3. * @param[in] b The 3-element vector representing the right-hand side of the equation. * @param[out] x The 3-element vector that will store the solution to the system. - * @return int that sepcifies whether the solve succeeded (1) or not (0). + * @return bool that sepcifies whether the solve succeeded (1) or not (0). */ template< typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE > GEOS_HOST_DEVICE inline -int solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) +bool solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE && x ) { LvArray::tensorOps::internal::checkSizes< 3, 3 >( A ); LvArray::tensorOps::internal::checkSizes< 3 >( b ); @@ -109,7 +109,7 @@ int solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE real64 const detA = LvArray::tensorOps::determinant< 3 >( A ); if( LvArray::math::abs( detA ) < singularMatrixTolerance ) - return 0; + return false; real64 const invA = 1.0 / detA; @@ -129,7 +129,7 @@ int solveThreeByThreeSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_TYPE x[1] = detX1 * invA; x[2] = detX2 * invA; - return 1; + return true; } /** @@ -169,17 +169,17 @@ void solveUpperTriangularSystem( MATRIX_TYPE const & A, RHS_TYPE const & b, SOL_ * @brief Solves a linear system using Gaussian elimination. * * This function performs Gaussian elimination on the given matrix `A` and right-hand side vector `b`. - * It transforms the matrix `A` into an upper triangular matrix and then solves for the solution `x` + * It transforms the matrix `A` boolo an upper triangular matrix and then solves for the solution `x` * using back substitution. * * @tparam N The size of the square matrix `A`. * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. * @tparam SOL_TYPE The type of the solution vector `x`. - * @param[in,out] A The matrix to be transformed into an upper triangular matrix. Modified in place. + * @param[in,out] A The matrix to be transformed boolo an upper triangular matrix. Modified in place. * @param[in,out] b The right-hand side vector. Modified in place to reflect the transformed system. * @param[out] x The solution vector. The result of solving the system `Ax = b`. - * @return int that sepcifies whether the solve succeeded (1) or not (0). + * @return bool that sepcifies whether the solve succeeded (1) or not (0). */ template< std::ptrdiff_t N, typename MATRIX_TYPE, @@ -187,14 +187,14 @@ template< std::ptrdiff_t N, typename SOL_TYPE > GEOS_HOST_DEVICE inline -int solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +bool solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); LvArray::tensorOps::internal::checkSizes< N, N >( A ); LvArray::tensorOps::internal::checkSizes< N >( b ); LvArray::tensorOps::internal::checkSizes< N >( x ); - // Step 1: Transform into an upper triangular matrix + // Step 1: Transform boolo an upper triangular matrix // 1.a. Find the pivot for( std::ptrdiff_t i = 0; i < N; ++i ) @@ -223,7 +223,7 @@ int solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) if( LvArray::math::abs( A[i][i] ) < singularMatrixTolerance ) - return 0; + return false; // 1.c Eliminate entries below the pivot for( std::ptrdiff_t k = i + 1; k < N; ++k ) @@ -240,7 +240,7 @@ int solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) // Step 2: Backward substitution solveUpperTriangularSystem< N >( A, b, std::forward< SOL_TYPE >( x ) ); - return 1; + return true; } } // details namespace @@ -256,22 +256,22 @@ int solveGaussianElimination( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) * @tparam MATRIX_TYPE The type of the matrix `A`. * @tparam RHS_TYPE The type of the right-hand side vector `b`. * @tparam SOL_TYPE The type of the solution vector `x`. - * @tparam MODIFY_MATRIX intean flag indicating whether the input matrix `A` and vector `b` should be modified. + * @tparam MODIFY_MATRIX boolean flag indicating whether the input matrix `A` and vector `b` should be modified. * If `1`, the matrix `A` and vector `b` are modified in place. If `0`, copies of * `A` and `b` are made, and the original data is left unchanged. * @param[in] A The matrix representing the coefficients of the system. * @param[in] b The right-hand side vector. * @param[out] x The solution vector. The result of solving the system `Ax = b`. - * @return int that sepcifies whether the solve succeeded (1) or not (0). + * @return bool that sepcifies whether the solve succeeded (1) or not (0). */ template< std::ptrdiff_t N, typename MATRIX_TYPE, typename RHS_TYPE, typename SOL_TYPE, - int MODIFY_MATRIX = 1 > + bool MODIFY_MATRIX = 1 > GEOS_HOST_DEVICE inline -int solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) +bool solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x ) { static_assert( N > 0, "N must be greater than 0." ); static_assert( N < 10, "N cannot be larger than 9" ); diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp index b0cd4c1eb4e..6f8bbd08bbd 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp @@ -48,21 +48,6 @@ class InvertibleLinearSystem : public LinearSystem< N > using LinearSystem< N >::solution; using LinearSystem< N >::rhs; - /// Default copy constructor - InvertibleLinearSystem( InvertibleLinearSystem const & ) = default; - - /// Default move constructor - InvertibleLinearSystem( InvertibleLinearSystem && ) = default; - - /// Deleted default constructor - InvertibleLinearSystem() = delete; - - /// Deleted copy assignment operator - InvertibleLinearSystem & operator=( InvertibleLinearSystem const & ) = delete; - - /// Deleted move assignment operator - InvertibleLinearSystem & operator=( InvertibleLinearSystem && ) = delete; - InvertibleLinearSystem( short int seed ) { std::mt19937 generator( seed ); @@ -99,21 +84,6 @@ class SingularLinearSystem : public LinearSystem< N > using LinearSystem< N >::solution; using LinearSystem< N >::rhs; - /// Default copy constructor - SingularLinearSystem( SingularLinearSystem const & ) = default; - - /// Default move constructor - SingularLinearSystem( SingularLinearSystem && ) = default; - - /// Deleted default constructor - SingularLinearSystem() = delete; - - /// Deleted copy assignment operator - SingularLinearSystem & operator=( SingularLinearSystem const & ) = delete; - - /// Deleted move assignment operator - SingularLinearSystem & operator=( SingularLinearSystem && ) = delete; - SingularLinearSystem( short int seed ) { std::mt19937 generator( seed ); @@ -167,7 +137,7 @@ class DenseLinearSolverTest : public ::testing::Test LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - int const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); PORTABLE_EXPECT_TRUE( success ); @@ -186,13 +156,13 @@ class DenseLinearSolverTest : public ::testing::Test forAll< parallelDevicePolicy<> >( 1, [=] GEOS_HOST_DEVICE ( int ) { - real64 sol[size]; - real64 matrix[size][size]; - real64 rhs[size]; + real64 sol[size]{}; + real64 matrix[size][size]{}; + real64 rhs[size]{}; LvArray::tensorOps::copy< size, size >( matrix, LS.matrix ); LvArray::tensorOps::copy< size >( rhs, LS.rhs ); - int const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); + bool const success = denseLinearAlgebra::solve< size >( matrix, rhs, sol ); PORTABLE_EXPECT_FALSE( success ); } ); @@ -232,7 +202,7 @@ TYPED_TEST_SUITE( DenseLinearSolverTest, Dimensions, NameGenerator ); TYPED_TEST( DenseLinearSolverTest, testDenseLA ) { this->test_solve(); - // this->test_singularSystem(); + this->test_singularSystem(); } int main( int argc, char * * argv ) diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp index 82a94ca86f3..17f093a569b 100644 --- a/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp +++ b/src/coreComponents/denseLinearAlgebra/unitTests/testUtils.hpp @@ -17,8 +17,8 @@ namespace testing #define PORTABLE_EXPECT_EQ( L, R ) GEOS_ERROR_IF_NE( L, R ) #define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) LVARRAY_ERROR_IF_GE_MSG( LvArray::math::abs( ( L ) -( R ) ), EPSILON, \ STRINGIZE( L ) " = " << ( L ) << "\n" << STRINGIZE( R ) " = " << ( R ) ); -#define PORTABLE_EXPECT_TRUE( value ) LVARRAY_ERROR_IF( value==0, "should be true" ) -#define PORTABLE_EXPECT_FALSE( value ) LVARRAY_ERROR_IF( value==1, "should be false" ) +#define PORTABLE_EXPECT_TRUE( value ) GEOS_ERROR_IF( !value, "should be true" ) +#define PORTABLE_EXPECT_FALSE( value ) GEOS_ERROR_IF( value, "should be false" ) #else #define PORTABLE_EXPECT_EQ( L, R ) EXPECT_EQ( L, R ) #define PORTABLE_EXPECT_NEAR( L, R, EPSILON ) EXPECT_LE( LvArray::math::abs( ( L ) -( R ) ), EPSILON ) << \