Skip to content

Commit

Permalink
fix bugs and add first two unitTests.
Browse files Browse the repository at this point in the history
  • Loading branch information
@CusiniM
CusiniM committed Aug 20, 2024
1 parent 026017a commit 640d56b
Show file tree
Hide file tree
Showing 3 changed files with 153 additions and 40 deletions.
82 changes: 43 additions & 39 deletions src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@

#include "common/DataTypes.hpp"
#include "denseLinearAlgebra/common/layouts.hpp"
#include "LvArray/src/tensorOps.hpp"
#include "common/logger/Logger.hpp"

#include <complex>

Expand All @@ -30,7 +32,7 @@ namespace geos
namespace denseLinearAlgebra
{

namespace internal
namespace details
{
/**
* @brief Solves a 2x2 linear system A * x = b.
Expand Down Expand Up @@ -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;
}

/**
Expand Down Expand Up @@ -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;
}

/**
Expand All @@ -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 )
{
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand All @@ -247,38 +251,38 @@ 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 )
{
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 ) );
}
}
}
Expand Down
3 changes: 2 additions & 1 deletion src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
set( serial_tests
testBlasLapack.cpp
testSolveLinearSystem.cpp )
testSolveLinearSystem.cpp
testDenseLASolvers.cpp )

set( dependencyList gtest denseLinearAlgebra )

Expand Down
108 changes: 108 additions & 0 deletions src/coreComponents/denseLinearAlgebra/unitTests/testDenseLASolvers.cpp
View file
Original file line number Diff line number Diff line change
@@ -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 <gtest/gtest.h>


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;
}

0 comments on commit 640d56b

Please sign in to comment.