diff --git a/include/micm/jit/solver/jit_lu_decomposition.inl b/include/micm/jit/solver/jit_lu_decomposition.inl index bb061030a..452e8c03f 100644 --- a/include/micm/jit/solver/jit_lu_decomposition.inl +++ b/include/micm/jit/solver/jit_lu_decomposition.inl @@ -89,7 +89,8 @@ namespace micm func.SetArrayElement(func.arguments_[2], U_ptr_index, JitType::Double, A_val); func.EndLoop(loop); } - else { + else + { auto loop = func.StartLoop("Uik_eq_zero_loop", 0, L); llvm::Value *zero_val = llvm::ConstantFP::get(*(func.context_), llvm::APFloat(0.0)); llvm::Value *iUf = llvm::ConstantInt::get(*(func.context_), llvm::APInt(64, uik_nkj->first)); @@ -146,7 +147,8 @@ namespace micm func.SetArrayElement(func.arguments_[1], L_ptr_index, JitType::Double, A_val); func.EndLoop(loop); } - else { + else + { auto loop = func.StartLoop("Lki_eq_zero_loop", 0, L); llvm::Value *zero_val = llvm::ConstantFP::get(*(func.context_), llvm::APFloat(0.0)); llvm::Value *iLf = llvm::ConstantInt::get(*(func.context_), llvm::APInt(64, lki_nkj->first)); diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl index 331eae0fe..6385e53e0 100644 --- a/include/micm/solver/lu_decomposition.inl +++ b/include/micm/solver/lu_decomposition.inl @@ -194,10 +194,12 @@ namespace micm // Upper trianglur matrix for (std::size_t iU = 0; iU < inLU.second; ++iU) { - if (*(do_aik++)){ + if (*(do_aik++)) + { U_vector[uik_nkj->first] = A_vector[*(aik++)]; } - else { + else + { U_vector[uik_nkj->first] = 0; } for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj) @@ -211,10 +213,12 @@ namespace micm L_vector[(lki_nkj++)->first] = 1.0; for (std::size_t iL = 0; iL < inLU.first; ++iL) { - if (*(do_aki++)){ + if (*(do_aki++)) + { L_vector[lki_nkj->first] = A_vector[*(aki++)]; } - else { + else + { L_vector[lki_nkj->first] = 0; } for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj) @@ -283,7 +287,8 @@ namespace micm std::copy(A_vector + *aik, A_vector + *aik + n_cells, U_vector + uik_nkj_first); ++aik; } - else { + else + { std::fill(U_vector + uik_nkj_first, U_vector + uik_nkj_first + n_cells, 0); } for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj) @@ -308,7 +313,8 @@ namespace micm std::copy(A_vector + *aki, A_vector + *aki + n_cells, L_vector + lki_nkj_first); ++aki; } - else { + else + { std::fill(L_vector + lki_nkj_first, L_vector + lki_nkj_first + n_cells, 0); } for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj) diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index fd9f850a1..3b4127a57 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -23,7 +23,7 @@ namespace micm const double h_max = parameters_.h_max_ == 0.0 ? time_step : std::min(time_step, parameters_.h_max_); const double h_start = parameters_.h_start_ == 0.0 ? std::max(parameters_.h_min_, DELTA_MIN) : std::min(h_max, parameters_.h_start_); - + SolverStats stats; double present_time = 0.0; diff --git a/include/micm/util/vector_matrix.hpp b/include/micm/util/vector_matrix.hpp index dbc9c1a9c..8fb7bc008 100644 --- a/include/micm/util/vector_matrix.hpp +++ b/include/micm/util/vector_matrix.hpp @@ -228,7 +228,8 @@ namespace micm return L; } - void Print() const { + void Print() const + { for (std::size_t i = 0; i < x_dim_; ++i) { for (std::size_t j = 0; j < y_dim_; ++j) diff --git a/src/solver/lu_decomposition.cu b/src/solver/lu_decomposition.cu index 90301c1e1..aef3b957d 100644 --- a/src/solver/lu_decomposition.cu +++ b/src/solver/lu_decomposition.cu @@ -62,7 +62,8 @@ namespace micm size_t A_idx = d_aik[aik_offset++] + tid; d_U[U_idx] = d_A[A_idx]; } - else { + else + { d_U[U_idx] = 0; } @@ -87,7 +88,8 @@ namespace micm size_t A_idx = d_aki[aki_offset++] + tid; d_L[L_idx] = d_A[A_idx]; } - else { + else + { d_L[L_idx] = 0; } diff --git a/test/unit/cuda/solver/test_cuda_linear_solver.cpp b/test/unit/cuda/solver/test_cuda_linear_solver.cpp index 81a4cf817..d2c714832 100644 --- a/test/unit/cuda/solver/test_cuda_linear_solver.cpp +++ b/test/unit/cuda/solver/test_cuda_linear_solver.cpp @@ -67,8 +67,7 @@ std::vector linearSolverGenerator(std::size_t number_of_blocks) CopyToDeviceDense(x); LinearSolverPolicy solver = LinearSolverPolicy(A, 0); - std::pair lu = - micm::LuDecomposition::GetLUMatrices(A, 0); + std::pair lu = micm::LuDecomposition::GetLUMatrices(A, 0); SparseMatrixPolicy lower_matrix = std::move(lu.first); SparseMatrixPolicy upper_matrix = std::move(lu.second); @@ -142,13 +141,21 @@ TEST(CudaLinearSolver, RandomMatrixVectorOrderingForGPU) TEST(CudaLinearSolver, AgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto initial_value : initial_values) + for (auto initial_value : initial_values) { - testExtremeInitialValue>(1, initial_value); - testExtremeInitialValue>(20, initial_value); - testExtremeInitialValue>( - 300, initial_value); - testExtremeInitialValue>( - 4000, initial_value); + testExtremeInitialValue>( + 1, initial_value); + testExtremeInitialValue< + Group20CudaDenseMatrix, + Group20CudaSparseMatrix, + micm::CudaLinearSolver>(20, initial_value); + testExtremeInitialValue< + Group300CudaDenseMatrix, + Group300CudaSparseMatrix, + micm::CudaLinearSolver>(300, initial_value); + testExtremeInitialValue< + Group4000CudaDenseMatrix, + Group4000CudaSparseMatrix, + micm::CudaLinearSolver>(4000, initial_value); } } \ No newline at end of file diff --git a/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp b/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp index a8a57d541..132cd30ab 100644 --- a/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp +++ b/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp @@ -29,13 +29,14 @@ void testCudaRandomMatrix(size_t n_grids) CPUSparseMatrixPolicy cpu_A(builder); GPUSparseMatrixPolicy gpu_A(builder); - // for nvhpc, the lognormal distribution produces significantly different values + // for nvhpc, the lognormal distribution produces significantly different values // for very large numbers of grid cells // To keep the accuracy on the check results function small, we only generat 1 blocks worth of // random values and then copy that into every other block for (std::size_t i = 0; i < 10; ++i) for (std::size_t j = 0; j < 10; ++j) - if (!cpu_A.IsZero(i, j)) { + if (!cpu_A.IsZero(i, j)) + { cpu_A[0][i][j] = get_double(); gpu_A[0][i][j] = cpu_A[0][i][j]; for (std::size_t i_block = 1; i_block < n_grids; ++i_block) @@ -43,7 +44,7 @@ void testCudaRandomMatrix(size_t n_grids) cpu_A[i_block][i][j] = cpu_A[0][i][j]; gpu_A[i_block][i][j] = cpu_A[0][i][j]; } - } + } micm::CudaLuDecomposition gpu_lud(gpu_A); auto gpu_LU = micm::CudaLuDecomposition::GetLUMatrices(gpu_A, 0); @@ -72,13 +73,13 @@ void testCudaRandomMatrix(size_t n_grids) { auto gpu_L = gpu_L_vector[i]; auto cpu_L = cpu_L_vector[i]; - EXPECT_LT(std::abs((gpu_L-cpu_L)/cpu_L), 1.0e-10); + EXPECT_LT(std::abs((gpu_L - cpu_L) / cpu_L), 1.0e-10); }; for (int j = 0; j < U_size; ++j) { auto gpu_U = gpu_U_vector[j]; auto cpu_U = cpu_U_vector[j]; - EXPECT_LT(std::abs((gpu_U-cpu_U)/cpu_U), 1.0e-10); + EXPECT_LT(std::abs((gpu_U - cpu_U) / cpu_U), 1.0e-10); }; } @@ -103,7 +104,8 @@ TEST(CudaLuDecomposition, RandomMatrixVectorOrdering) TEST(CudaLuDecomposition, AgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto& value : initial_values) { + for (auto& value : initial_values) + { testExtremeValueInitialization(1, value); testExtremeValueInitialization(100, value); testExtremeValueInitialization(1000, value); diff --git a/test/unit/jit/solver/test_jit_linear_solver.cpp b/test/unit/jit/solver/test_jit_linear_solver.cpp index 795b73152..0ea39df87 100644 --- a/test/unit/jit/solver/test_jit_linear_solver.cpp +++ b/test/unit/jit/solver/test_jit_linear_solver.cpp @@ -47,7 +47,11 @@ TEST(JitLinearSolver, DiagonalMatrixVectorOrdering) TEST(JitLinearSolver, AgnosticToInitialValue) { double initial_values[5] = { -INFINITY, INFINITY }; - for(auto initial_value : initial_values) { - testExtremeInitialValue>(1, initial_value); + for (auto initial_value : initial_values) + { + testExtremeInitialValue< + Group1VectorMatrix, + Group1SparseVectorMatrix, + micm::JitLinearSolver<1, Group1SparseVectorMatrix>>(1, initial_value); } } \ No newline at end of file diff --git a/test/unit/jit/solver/test_jit_lu_decomposition.cpp b/test/unit/jit/solver/test_jit_lu_decomposition.cpp index 7901f4855..531c04ab2 100644 --- a/test/unit/jit/solver/test_jit_lu_decomposition.cpp +++ b/test/unit/jit/solver/test_jit_lu_decomposition.cpp @@ -46,7 +46,8 @@ TEST(JitLuDecomposition, DiagonalMatrixVectorOrdering) TEST(JitLuDecomposition, AgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto& value : initial_values) { + for (auto& value : initial_values) + { testExtremeValueInitialization>(1, value); testExtremeValueInitialization>(2, value); testExtremeValueInitialization>(3, value); diff --git a/test/unit/solver/test_linear_solver.cpp b/test/unit/solver/test_linear_solver.cpp index 4062220b9..bf1e7775f 100644 --- a/test/unit/solver/test_linear_solver.cpp +++ b/test/unit/solver/test_linear_solver.cpp @@ -38,7 +38,7 @@ TEST(LinearSolver, DiagonalMarkowitzReorder) TEST(LinearSolver, StandardOrderingAgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto initial_value : initial_values) + for (auto initial_value : initial_values) testExtremeInitialValue>(5, initial_value); } @@ -71,11 +71,16 @@ TEST(LinearSolver, RandomMatrixVectorOrdering) TEST(LinearSolver, VectorOrderingAgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto initial_value : initial_values) { - testExtremeInitialValue>(1, initial_value); - testExtremeInitialValue>(2, initial_value); - testExtremeInitialValue>(5, initial_value); - testExtremeInitialValue>(5, initial_value); + for (auto initial_value : initial_values) + { + testExtremeInitialValue>( + 1, initial_value); + testExtremeInitialValue>( + 2, initial_value); + testExtremeInitialValue>( + 5, initial_value); + testExtremeInitialValue>( + 5, initial_value); } } diff --git a/test/unit/solver/test_linear_solver_policy.hpp b/test/unit/solver/test_linear_solver_policy.hpp index a8db34e16..56b65131c 100644 --- a/test/unit/solver/test_linear_solver_policy.hpp +++ b/test/unit/solver/test_linear_solver_policy.hpp @@ -219,9 +219,12 @@ void testExtremeInitialValue(std::size_t number_of_blocks, double initial_value) const size_t size = 30; auto builder = SparseMatrixPolicy::Create(size).SetNumberOfBlocks(number_of_blocks).InitialValue(0); - for (std::size_t i = 0; i < size; ++i) { - for (std::size_t j = 0; j < size; ++j) { - if (i == j || gen_bool()) { + for (std::size_t i = 0; i < size; ++i) + { + for (std::size_t j = 0; j < size; ++j) + { + if (i == j || gen_bool()) + { builder = builder.WithElement(i, j); } } diff --git a/test/unit/solver/test_lu_decomposition.cpp b/test/unit/solver/test_lu_decomposition.cpp index fff36d369..13b089b47 100644 --- a/test/unit/solver/test_lu_decomposition.cpp +++ b/test/unit/solver/test_lu_decomposition.cpp @@ -37,7 +37,8 @@ TEST(LuDecomposition, DiagonalMatrixStandardOrdering) TEST(LuDecomposition, AgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto& value : initial_values) { + for (auto& value : initial_values) + { testExtremeValueInitialization(5, value); } } @@ -77,7 +78,8 @@ TEST(LuDecomposition, DiagonalMatrixVectorOrdering) TEST(LuDecomposition, VectorOrderingAgnosticToInitialValue) { double initial_values[5] = { -INFINITY, -1.0, 0.0, 1.0, INFINITY }; - for(auto& value : initial_values) { + for (auto& value : initial_values) + { testExtremeValueInitialization(5, value); testExtremeValueInitialization(5, value); testExtremeValueInitialization(5, value); diff --git a/test/unit/solver/test_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp index b62838652..9c9067a00 100644 --- a/test/unit/solver/test_lu_decomposition_policy.hpp +++ b/test/unit/solver/test_lu_decomposition_policy.hpp @@ -132,8 +132,7 @@ template void testSingularMatrix() { SparseMatrixPolicy A = SparseMatrixPolicy( - SparseMatrixPolicy::Create(2).InitialValue(0).WithElement(0, 0).WithElement(0, 1).WithElement(1, 0).WithElement( - 1, 1)); + SparseMatrixPolicy::Create(2).InitialValue(0).WithElement(0, 0).WithElement(0, 1).WithElement(1, 0).WithElement(1, 1)); A[0][0][0] = 0; A[0][0][1] = 1; @@ -175,9 +174,7 @@ void testRandomMatrix(std::size_t number_of_blocks) bool is_singular{ false }; lud.template Decompose(A, LU.first, LU.second, is_singular); check_results( - A, LU.first, LU.second, [&](const double a, const double b) -> void { - EXPECT_NEAR(a, b, 1.0e-9); - }); + A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-9); }); } template @@ -195,13 +192,14 @@ void testExtremeValueInitialization(std::size_t number_of_blocks, double initial SparseMatrixPolicy A(builder); - // for nvhpc, the lognormal distribution produces significantly different values + // for nvhpc, the lognormal distribution produces significantly different values // for very large numbers of grid cells // To keep the accuracy on the check results function small, we only generat 1 blocks worth of // random values and then copy that into every other block for (std::size_t i = 0; i < size; ++i) for (std::size_t j = 0; j < size; ++j) - if (!A.IsZero(i, j)){ + if (!A.IsZero(i, j)) + { A[0][i][j] = get_double(); for (std::size_t i_block = 1; i_block < number_of_blocks; ++i_block) A[i_block][i][j] = A[0][i][j]; @@ -223,9 +221,7 @@ void testExtremeValueInitialization(std::size_t number_of_blocks, double initial CopyToHost(LU.second); check_results( - A, LU.first, LU.second, [&](const double a, const double b) -> void { - EXPECT_NEAR(a, b, 1.0e-09); - }); + A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-09); }); } template