From d8bdb034e5a797adea82f2e375a2abae512ed048 Mon Sep 17 00:00:00 2001 From: jsandham Date: Sat, 25 Jul 2026 13:57:44 -0400 Subject: [PATCH 1/3] Add ruiz scaling to csr_matrix --- clients/testing/CMakeLists.txt | 1 + clients/testing/test_arguments.h | 2 +- clients/testing/test_enums.h | 15 ++- clients/testing/test_functions.h | 3 +- clients/testing/test_functions_dispatch.cpp | 6 +- .../testing/test_functions_ruiz_scaling.cpp | 89 +++++++++++++ clients/testing/test_yaml_loader.h | 4 +- clients/testing/tests/CMakeLists.txt | 1 + clients/testing/tests/test_exclusive_scan.cpp | 4 +- clients/testing/tests/test_ruiz_scaling.cpp | 29 +++++ clients/testing/tests/test_ruiz_scaling.yaml | 29 +++++ .../testing/tests/test_tridiagonal_solver.cpp | 2 +- library/include/csr_matrix.h | 2 + library/src/backend/device/CMakeLists.txt | 1 + .../src/backend/device/cuda/CMakeLists.txt | 1 + library/src/backend/device/cuda/cuda_math.h | 1 + .../backend/device/cuda/cuda_ruiz_scaling.cu | 46 +++++++ .../backend/device/cuda/cuda_ruiz_scaling.h | 36 ++++++ library/src/backend/device/device_math.h | 1 + .../backend/device/device_ruiz_scaling.cpp | 51 ++++++++ .../src/backend/device/device_ruiz_scaling.h | 37 ++++++ library/src/backend/host/CMakeLists.txt | 1 + library/src/backend/host/host_math.h | 1 + .../src/backend/host/host_ruiz_scaling.cpp | 118 ++++++++++++++++++ library/src/backend/host/host_ruiz_scaling.h | 38 ++++++ library/src/csr_matrix.cpp | 12 ++ 26 files changed, 516 insertions(+), 15 deletions(-) create mode 100644 clients/testing/test_functions_ruiz_scaling.cpp create mode 100644 clients/testing/tests/test_ruiz_scaling.cpp create mode 100644 clients/testing/tests/test_ruiz_scaling.yaml create mode 100644 library/src/backend/device/cuda/cuda_ruiz_scaling.cu create mode 100644 library/src/backend/device/cuda/cuda_ruiz_scaling.h create mode 100644 library/src/backend/device/device_ruiz_scaling.cpp create mode 100644 library/src/backend/device/device_ruiz_scaling.h create mode 100644 library/src/backend/host/host_ruiz_scaling.cpp create mode 100644 library/src/backend/host/host_ruiz_scaling.h diff --git a/clients/testing/CMakeLists.txt b/clients/testing/CMakeLists.txt index 129bc44..9f75163 100644 --- a/clients/testing/CMakeLists.txt +++ b/clients/testing/CMakeLists.txt @@ -92,6 +92,7 @@ target_sources(test_main test_functions_spgeam.cpp test_functions_transpose.cpp test_functions_transpose_dense.cpp + test_functions_ruiz_scaling.cpp test_functions_csric0.cpp test_functions_csrilu0.cpp test_functions_tridiagonal_solver.cpp diff --git a/clients/testing/test_arguments.h b/clients/testing/test_arguments.h index 520b519..5c75efb 100644 --- a/clients/testing/test_arguments.h +++ b/clients/testing/test_arguments.h @@ -38,7 +38,7 @@ namespace testing struct Arguments { testing::category category; // iterative_solvers, math, primitive - testing::fixture fixture; // Jacobi, CG, SpMV, ExclusiveScan, etc + testing::fixture fixture; // Jacobi, CG, SpMV, exclusive_scan, etc std::string group; // small, medium, large, etc std::string filename; // bmwcra_1.mtx, shipsec1.mtx, etc testing::backend backend; // CPU, GPU diff --git a/clients/testing/test_enums.h b/clients/testing/test_enums.h index 4c6302b..62af2ab 100644 --- a/clients/testing/test_enums.h +++ b/clients/testing/test_enums.h @@ -67,8 +67,9 @@ namespace testing transpose_dense, CSRIC0, CSRILU0, - TridiagonalSolver, - ExclusiveScan, + tridiagonal_solver, + exclusive_scan, + ruiz_scaling, Unknown }; @@ -212,10 +213,12 @@ namespace testing return "CSRIC0"; case fixture::CSRILU0: return "CSRILU0"; - case fixture::TridiagonalSolver: - return "TridiagonalSolver"; - case fixture::ExclusiveScan: - return "ExclusiveScan"; + case fixture::tridiagonal_solver: + return "tridiagonal_solver"; + case fixture::exclusive_scan: + return "exclusive_scan"; + case fixture::ruiz_scaling: + return "ruiz_scaling"; } return "Invalid"; diff --git a/clients/testing/test_functions.h b/clients/testing/test_functions.h index f59f2dc..1561332 100644 --- a/clients/testing/test_functions.h +++ b/clients/testing/test_functions.h @@ -2,7 +2,7 @@ // // MIT License // -// Copyright(c) 2025 James Sandham +// Copyright(c) 2025-2026 James Sandham // // Permission is hereby granted, free of charge, to any person obtaining a copy // of this softwareand associated documentation files(the "Software"), to deal @@ -51,6 +51,7 @@ namespace testing bool test_exclusive_scan(Arguments arg); // csr matrix + bool test_ruiz_scaling(Arguments arg); bool test_transpose(Arguments arg); bool test_transpose_dense(Arguments arg); bool test_multiply_by_vector(Arguments arg); diff --git a/clients/testing/test_functions_dispatch.cpp b/clients/testing/test_functions_dispatch.cpp index 6630567..a3b141f 100644 --- a/clients/testing/test_functions_dispatch.cpp +++ b/clients/testing/test_functions_dispatch.cpp @@ -72,7 +72,7 @@ namespace testing { switch(arg.fixture) { - case fixture::TridiagonalSolver: + case fixture::tridiagonal_solver: return test_tridiagonal_solver(arg); } @@ -100,7 +100,7 @@ namespace testing { switch(arg.fixture) { - case fixture::ExclusiveScan: + case fixture::exclusive_scan: return test_exclusive_scan(arg); } @@ -127,6 +127,8 @@ namespace testing return test_transpose(arg); case fixture::transpose_dense: return test_transpose_dense(arg); + case fixture::ruiz_scaling: + return test_ruiz_scaling(arg); } return false; } diff --git a/clients/testing/test_functions_ruiz_scaling.cpp b/clients/testing/test_functions_ruiz_scaling.cpp new file mode 100644 index 0000000..5cd7084 --- /dev/null +++ b/clients/testing/test_functions_ruiz_scaling.cpp @@ -0,0 +1,89 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include "test_functions.h" +#include "utility.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "linalg.h" + +bool testing::test_ruiz_scaling(Arguments arg) +{ + linalg::csr_matrix mat_A; + mat_A.read_mtx(arg.filename); + + linalg::csr_matrix mat_A_copy; + mat_A_copy.copy_from(mat_A); + + linalg::vector D1(mat_A.get_m()); + linalg::vector D2(mat_A.get_m()); + + if(arg.backend == backend::GPU) + { + mat_A.move_to_device(); + mat_A_copy.move_to_device(); + D1.move_to_device(); + D2.move_to_device(); + } + + // Warmup + for(int i = 0; i < 4; i++) + { + mat_A.apply_ruiz_scaling(D1, D2); + mat_A.copy_from(mat_A_copy); + } + linalg::synchronize(); + + // Timed run + auto t1 = std::chrono::high_resolution_clock::now(); + for(int i = 0; i < 100; i++) + { + mat_A.apply_ruiz_scaling(D1, D2); + mat_A.copy_from(mat_A_copy); + } + linalg::synchronize(); + auto t2 = std::chrono::high_resolution_clock::now(); + + std::chrono::duration ms_float = t2 - t1; + std::cout << "Solve time: " << ms_float.count() << "ms" << std::endl; + + if(arg.backend == backend::GPU) + { + mat_A.move_to_host(); + D1.move_to_host(); + D2.move_to_host(); + } + + return true; +} diff --git a/clients/testing/test_yaml_loader.h b/clients/testing/test_yaml_loader.h index 58d0b0b..f6e3233 100644 --- a/clients/testing/test_yaml_loader.h +++ b/clients/testing/test_yaml_loader.h @@ -404,8 +404,8 @@ inline testing::fixture string_to_fixture(const std::string& str) {"transpose_dense", testing::fixture::transpose_dense}, {"CSRIC0", testing::fixture::CSRIC0}, {"CSRILU0", testing::fixture::CSRILU0}, - {"TridiagonalSolver", testing::fixture::TridiagonalSolver}, - {"ExclusiveScan", testing::fixture::ExclusiveScan}}; + {"tridiagonal_solver", testing::fixture::tridiagonal_solver}, + {"exclusive_scan", testing::fixture::exclusive_scan}}; // Find the string in the map auto it = fixture_map.find(str); diff --git a/clients/testing/tests/CMakeLists.txt b/clients/testing/tests/CMakeLists.txt index 47f0b1e..a918719 100644 --- a/clients/testing/tests/CMakeLists.txt +++ b/clients/testing/tests/CMakeLists.txt @@ -49,4 +49,5 @@ target_sources(test_main test_transpose_dense.cpp test_tridiagonal_solver.cpp test_exclusive_scan.cpp + test_ruiz_scaling.cpp ) diff --git a/clients/testing/tests/test_exclusive_scan.cpp b/clients/testing/tests/test_exclusive_scan.cpp index d056ca6..1316cd1 100644 --- a/clients/testing/tests/test_exclusive_scan.cpp +++ b/clients/testing/tests/test_exclusive_scan.cpp @@ -2,7 +2,7 @@ // // MIT License // -// Copyright(c) 2025 James Sandham +// Copyright(c) 2025-2026 James Sandham // // Permission is hereby granted, free of charge, to any person obtaining a copy // of this softwareand associated documentation files(the "Software"), to deal @@ -26,4 +26,4 @@ #include "../test.h" -INSTANTIATE_TEST(primitive, ExclusiveScan, exclusive_scan, "tests/test_exclusive_scan.yaml"); +INSTANTIATE_TEST(primitive, exclusive_scan, exclusive_scan, "tests/test_exclusive_scan.yaml"); diff --git a/clients/testing/tests/test_ruiz_scaling.cpp b/clients/testing/tests/test_ruiz_scaling.cpp new file mode 100644 index 0000000..5c1e4d2 --- /dev/null +++ b/clients/testing/tests/test_ruiz_scaling.cpp @@ -0,0 +1,29 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include "../test.h" + +INSTANTIATE_TEST(csr_matrix, ruiz_scaling, ruiz_scaling, "tests/test_ruiz_scaling.yaml"); diff --git a/clients/testing/tests/test_ruiz_scaling.yaml b/clients/testing/tests/test_ruiz_scaling.yaml new file mode 100644 index 0000000..c36d9d0 --- /dev/null +++ b/clients/testing/tests/test_ruiz_scaling.yaml @@ -0,0 +1,29 @@ +Tests: + small: + matrix_file: ["matrices/SPD/bcsstm02/bcsstm02.mtx", + "matrices/SPD/ex5/ex5.mtx", + "matrices/SPD/mesh1em6/mesh1em6.mtx", + "matrices/SPD/bcsstk01/bcsstk01.mtx", + "matrices/SPD/mesh3em5/mesh3em5.mtx", + "matrices/SPD/mesh2em5/mesh2em5.mtx", + "matrices/SPD/fv1/fv1.mtx", + "matrices/SPD/fv2/fv2.mtx", + "matrices/SPD/fv3/fv3.mtx", + "matrices/SPD/bodyy4/bodyy4.mtx", + "matrices/SPD/shallow_water1/shallow_water1.mtx", + "matrices/SPD/shallow_water2/shallow_water2.mtx", + "matrices/SPD/gyro_m/gyro_m.mtx", + "matrices/SPD/crystm02/crystm02.mtx", + "matrices/SPD/crystm03/crystm03.mtx", + "matrices/SPD/finan512/finan512.mtx", + "matrices/SPD/thermomech_TC/thermomech_TC.mtx", + "matrices/SPD/thermomech_TK/thermomech_TK.mtx", + "matrices/SPD/bundle1/bundle1.mtx", + "matrices/SPD/nos6/nos6.mtx", + "matrices/SPD/bodyy6/bodyy6.mtx", + "matrices/SPD/cvxbqp1/cvxbqp1.mtx", + "matrices/SPD/bcsstk38/bcsstk38.mtx", + "matrices/SPD/wathen100/wathen100.mtx", + "matrices/SPD/apache1/apache1.mtx", + "matrices/SPD/wathen120/wathen120.mtx"] + backend: [CPU] diff --git a/clients/testing/tests/test_tridiagonal_solver.cpp b/clients/testing/tests/test_tridiagonal_solver.cpp index 46797cd..e4e5973 100644 --- a/clients/testing/tests/test_tridiagonal_solver.cpp +++ b/clients/testing/tests/test_tridiagonal_solver.cpp @@ -27,6 +27,6 @@ #include "../test.h" INSTANTIATE_TEST(direct_solvers, - TridiagonalSolver, + tridiagonal_solver, tridiagonal_solver, "tests/test_tridiagonal_solver.yaml"); diff --git a/library/include/csr_matrix.h b/library/include/csr_matrix.h index b4522a4..e8eda1a 100644 --- a/library/include/csr_matrix.h +++ b/library/include/csr_matrix.h @@ -365,6 +365,8 @@ namespace linalg */ void transpose(csr_matrix& T) const; + void apply_ruiz_scaling(vector& D1, vector& D2); + /*! \brief Reads a sparse matrix from a Matrix Market (.mtx) file into this CSR object. * \param filename The path to the .mtx file. * \return `true` if the file was successfully read and parsed, `false` otherwise. diff --git a/library/src/backend/device/CMakeLists.txt b/library/src/backend/device/CMakeLists.txt index b9d033c..dafcf84 100644 --- a/library/src/backend/device/CMakeLists.txt +++ b/library/src/backend/device/CMakeLists.txt @@ -34,6 +34,7 @@ target_sources(linalglib device_matrix_vector.cpp device_extract.cpp device_scale.cpp + device_ruiz_scaling.cpp device_ssor.cpp device_csr2csc.cpp device_csrtrsv.cpp diff --git a/library/src/backend/device/cuda/CMakeLists.txt b/library/src/backend/device/cuda/CMakeLists.txt index b6ea629..66586dd 100644 --- a/library/src/backend/device/cuda/CMakeLists.txt +++ b/library/src/backend/device/cuda/CMakeLists.txt @@ -33,6 +33,7 @@ target_sources(linalglib cuda_ssor.cu cuda_extract.cu cuda_scale.cu + cuda_ruiz_scaling.cu cuda_csric0.cu cuda_csrilu0.cu cuda_csr2csc.cu diff --git a/library/src/backend/device/cuda/cuda_math.h b/library/src/backend/device/cuda/cuda_math.h index 03b4fbf..2b4d5a2 100644 --- a/library/src/backend/device/cuda/cuda_math.h +++ b/library/src/backend/device/cuda/cuda_math.h @@ -35,6 +35,7 @@ #include "cuda_csrtrsv.h" #include "cuda_extract.h" #include "cuda_matrix_vector.h" +#include "cuda_ruiz_scaling.h" #include "cuda_scale.h" #include "cuda_ssor.h" #include "cuda_tridiagonal.h" diff --git a/library/src/backend/device/cuda/cuda_ruiz_scaling.cu b/library/src/backend/device/cuda/cuda_ruiz_scaling.cu new file mode 100644 index 0000000..06adb45 --- /dev/null +++ b/library/src/backend/device/cuda/cuda_ruiz_scaling.cu @@ -0,0 +1,46 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include + +#include "cuda_ruiz_scaling.h" + +#include "../../../trace.h" + +//------------------------------------------------------------------------------- +// Ruiz scaling +//------------------------------------------------------------------------------- +template +void linalg::cuda_ruiz_scaling(T* D1, + const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, T* D2) +{ + ROUTINE_TRACE("linalg::cuda_ruiz_scaling_impl"); + //scale_diagonal_kernel<256> + // <<<((m - 1) / 256 + 1), 256>>>(m, csr_row_ptr, csr_col_ind, csr_val, scalar); + //CHECK_CUDA_LAUNCH_ERROR(); +} + +template void linalg::cuda_ruiz_scaling(double* D1, const int* csr_row_ptr, const int* csr_col_ind, double* csr_val, int m, double* D2); diff --git a/library/src/backend/device/cuda/cuda_ruiz_scaling.h b/library/src/backend/device/cuda/cuda_ruiz_scaling.h new file mode 100644 index 0000000..8a447ae --- /dev/null +++ b/library/src/backend/device/cuda/cuda_ruiz_scaling.h @@ -0,0 +1,36 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** +#ifndef CUDA_RUIZ_sCALING_H +#define CUDA_RUIZ_SCALING_H + +namespace linalg +{ + template + void cuda_ruiz_scaling( + T* D1, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, T* D2); +} + +#endif diff --git a/library/src/backend/device/device_math.h b/library/src/backend/device/device_math.h index ebf4b1a..cc2a03c 100644 --- a/library/src/backend/device/device_math.h +++ b/library/src/backend/device/device_math.h @@ -38,6 +38,7 @@ #include "device_csrtrsv.h" #include "device_extract.h" #include "device_matrix_vector.h" +#include "device_ruiz_scaling.h" #include "device_scale.h" #include "device_ssor.h" #include "device_tridiagonal.h" diff --git a/library/src/backend/device/device_ruiz_scaling.cpp b/library/src/backend/device/device_ruiz_scaling.cpp new file mode 100644 index 0000000..d028179 --- /dev/null +++ b/library/src/backend/device/device_ruiz_scaling.cpp @@ -0,0 +1,51 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include "device_ruiz_scaling.h" + +#include + +#include "../../trace.h" +#include "../../utility.h" + +#if defined(LINALGLIB_HAS_CUDA) +#include "cuda/cuda_ruiz_scaling.h" +#endif + +void linalg::device_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2) +{ + ROUTINE_TRACE("linalg::device_ruiz_scaling"); + if constexpr(is_cuda_available()) + { + CALL_CUDA(cuda_ruiz_scaling( + D1.get_vec(), A.get_row_ptr(), A.get_col_ind(), A.get_val(), A.get_m(), D2.get_vec())); + } + else + { + std::cout << "Error: Not device backend available for the function " << __func__ + << std::endl; + } +} diff --git a/library/src/backend/device/device_ruiz_scaling.h b/library/src/backend/device/device_ruiz_scaling.h new file mode 100644 index 0000000..1d1381a --- /dev/null +++ b/library/src/backend/device/device_ruiz_scaling.h @@ -0,0 +1,37 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** +#ifndef DEVICE_RUIZ_SCALING_H +#define DEVICE_RUIZ_SCALING_H + +#include "csr_matrix.h" +#include "vector.h" + +namespace linalg +{ + void device_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2); +} + +#endif diff --git a/library/src/backend/host/CMakeLists.txt b/library/src/backend/host/CMakeLists.txt index 1320c75..f98b532 100644 --- a/library/src/backend/host/CMakeLists.txt +++ b/library/src/backend/host/CMakeLists.txt @@ -34,6 +34,7 @@ target_sources(linalglib host_matrix_vector.cpp host_extract.cpp host_scale.cpp + host_ruiz_scaling.cpp host_ssor.cpp host_csr2csc.cpp host_csrtrsv.cpp diff --git a/library/src/backend/host/host_math.h b/library/src/backend/host/host_math.h index d7494ae..302984f 100644 --- a/library/src/backend/host/host_math.h +++ b/library/src/backend/host/host_math.h @@ -38,6 +38,7 @@ #include "host_csrtrsv.h" #include "host_extract.h" #include "host_matrix_vector.h" +#include "host_ruiz_scaling.h" #include "host_scale.h" #include "host_ssor.h" #include "host_tridiagonal.h" diff --git a/library/src/backend/host/host_ruiz_scaling.cpp b/library/src/backend/host/host_ruiz_scaling.cpp new file mode 100644 index 0000000..3d21d94 --- /dev/null +++ b/library/src/backend/host/host_ruiz_scaling.cpp @@ -0,0 +1,118 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include +#include +#include + +#include "host_ruiz_scaling.h" + +#include "../../trace.h" + +namespace linalg +{ + template + static void host_ruiz_scaling_impl( + T* D1, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, T* D2) + { + ROUTINE_TRACE("host_scale_diagonal_impl"); + + int k = 10; + + // D_1^0 = I, D_2^0 = I + for(int i = 0; i < m; i++) + { + D1[i] = static_cast(1); + D2[i] = static_cast(1); + } + + std::vector DR(m, static_cast(0)); + std::vector DC(m, static_cast(0)); + + for(int iter = 0; iter < k; iter++) + { + for(int i = 0; i < m; i++) + { + const int start = csr_row_ptr[i]; + const int end = csr_row_ptr[i + 1]; + + for(int j = start; j < end; j++) + { + const int col = csr_col_ind[j]; + const T val = csr_val[j]; + + DR[i] = std::max(DR[i], std::abs(val)); + DC[col] = std::max(DC[col], std::abs(val)); + } + } + + T row_max = static_cast(0); + T col_max = static_cast(0); + for(int i = 0; i < m; i++) + { + row_max = std::max(row_max, DR[i]); + col_max = std::max(col_max, DC[i]); + + DR[i] = std::sqrt(DR[i]); + DC[i] = std::sqrt(DC[i]); + } + + std::cout << "row_max: " << row_max << " col_max: " << col_max << std::endl; + + // Update A_k+1 = DR^-1 * A_k * DC^-1 + for(int i = 0; i < m; i++) + { + const int start = csr_row_ptr[i]; + const int end = csr_row_ptr[i + 1]; + + for(int j = start; j < end; j++) + { + const int col = csr_col_ind[j]; + csr_val[j] /= (DR[i] * DC[col]); + } + } + + // Update D1_k+1 = D1_k * DR^-1, D2_k+1 = D2_k * DC^-1 + for(int i = 0; i < m; i++) + { + D1[i] = D1[i] / DR[i]; + D2[i] = D2[i] / DC[i]; + } + } + } +} + +void linalg::host_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2) +{ + ROUTINE_TRACE("linalg::host_ruiz_scaling"); + + assert(A.get_m() == A.get_n()); + assert(D1.get_size() == A.get_m()); + assert(D2.get_size() == A.get_m()); + + host_ruiz_scaling_impl( + D1.get_vec(), A.get_row_ptr(), A.get_col_ind(), A.get_val(), A.get_m(), D2.get_vec()); +} diff --git a/library/src/backend/host/host_ruiz_scaling.h b/library/src/backend/host/host_ruiz_scaling.h new file mode 100644 index 0000000..5c338a5 --- /dev/null +++ b/library/src/backend/host/host_ruiz_scaling.h @@ -0,0 +1,38 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#ifndef HOST_RUIZ_SCALING_H +#define HOST_RUIZ_SCALING_H + +#include "csr_matrix.h" +#include "vector.h" + +namespace linalg +{ + void host_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2); +} + +#endif // HOST_APPLY_RUIZ_SCALING_H diff --git a/library/src/csr_matrix.cpp b/library/src/csr_matrix.cpp index 9873611..6f1aa97 100644 --- a/library/src/csr_matrix.cpp +++ b/library/src/csr_matrix.cpp @@ -418,6 +418,18 @@ void csr_matrix::transpose(csr_matrix& T) const "linalg::csr_matrix::transpose", host_transpose_matrix, device_transpose_matrix, *this, T); } +void csr_matrix::apply_ruiz_scaling(vector& D1, vector& D2) +{ + ROUTINE_TRACE("csr_matrix::apply_ruiz_scaling"); + + backend_dispatch("linalg::csr_matrix::apply_ruiz_scaling", + host_ruiz_scaling, + device_ruiz_scaling, + D1, + *this, + D2); +} + // Structure to hold triplet (COO) format data struct triplet { From fda3c598cf6f0aa760c5bf9da09a84d4157e4781 Mon Sep 17 00:00:00 2001 From: jsandham Date: Wed, 29 Jul 2026 19:13:05 -0400 Subject: [PATCH 2/3] Add GPU ruiz symmetric scaling --- clients/examples/rsamg_example.cpp | 8 +- clients/examples/saamg_example.cpp | 9 +- clients/examples/uaamg_example.cpp | 5 +- clients/testing/CMakeLists.txt | 1 + clients/testing/test_arguments.h | 14 +- clients/testing/test_enums.h | 161 +++++++-------- clients/testing/test_functions.h | 1 + clients/testing/test_functions_amg.cpp | 32 +-- clients/testing/test_functions_classical.cpp | 12 +- clients/testing/test_functions_dispatch.cpp | 14 +- clients/testing/test_functions_krylov.cpp | 6 +- .../testing/test_functions_ruiz_scaling.cpp | 46 +++-- .../test_functions_symmetric_ruiz_scaling.cpp | 96 +++++++++ .../test_functions_triangular_solve.cpp | 12 +- .../test_functions_tridiagonal_solver.cpp | 4 +- clients/testing/test_yaml_loader.h | 134 ++++++------- clients/testing/tests/CMakeLists.txt | 1 + clients/testing/tests/test_BICGSTAB.yaml | 2 +- clients/testing/tests/test_CG.yaml | 16 +- clients/testing/tests/test_GMRES.yaml | 2 +- clients/testing/tests/test_gauss_seidel.cpp | 2 +- clients/testing/tests/test_jacobi.cpp | 2 +- clients/testing/tests/test_ruiz_scaling.yaml | 49 +++++ .../tests/test_symmetric_gauss_seidel.cpp | 2 +- .../tests/test_symmetric_ruiz_scaling.cpp | 32 +++ .../tests/test_symmetric_ruiz_scaling.yaml | 78 ++++++++ .../testing/tests/test_triangular_solve.yaml | 4 +- .../tests/test_tridiagonal_solver.yaml | 14 +- library/include/csr_matrix.h | 3 +- library/include/iterative_solvers/amg/amg.h | 38 ++-- .../backend/device/cuda/cuda_ruiz_scaling.cu | 92 ++++++++- .../backend/device/cuda/cuda_ruiz_scaling.h | 14 +- .../device/cuda/ruiz_scaling_kernels.cuh | 185 ++++++++++++++++++ .../backend/device/device_ruiz_scaling.cpp | 28 ++- .../src/backend/device/device_ruiz_scaling.h | 5 +- .../src/backend/host/host_ruiz_scaling.cpp | 172 ++++++++++++++-- library/src/backend/host/host_ruiz_scaling.h | 5 +- library/src/csr_matrix.cpp | 19 +- library/src/iterative_solvers/amg/amg.cpp | 30 +-- 39 files changed, 1040 insertions(+), 310 deletions(-) create mode 100644 clients/testing/test_functions_symmetric_ruiz_scaling.cpp create mode 100644 clients/testing/tests/test_symmetric_ruiz_scaling.cpp create mode 100644 clients/testing/tests/test_symmetric_ruiz_scaling.yaml create mode 100644 library/src/backend/device/cuda/ruiz_scaling_kernels.cuh diff --git a/clients/examples/rsamg_example.cpp b/clients/examples/rsamg_example.cpp index afa8b17..a91af00 100644 --- a/clients/examples/rsamg_example.cpp +++ b/clients/examples/rsamg_example.cpp @@ -50,9 +50,9 @@ int main() control.rel_tol = 1e-08; control.abs_tol = 1e-08; - // int cycles = linalg::amg_solve(hierachy, x.data(), b.data(), 10, 10, Cycle::Vcycle, Smoother::Gauss_Seidel, control); - // int cycles = linalg::amg_solve(hierachy, x.data(), b.data(), 2, 2, Cycle::Wcycle, Smoother::Gauss_Seidel, control); - // int cycles = linalg::amg_solve(hierachy, x.data(), b.data(), 2, 2, Cycle::Wcycle, Smoother::SOR, control); + // int cycles = linalg::amg_solve(hierachy, x.data(), b.data(), 10, 10, cycle::vcycle, smoother::gauss_seidel, control); + // int cycles = linalg::amg_solve(hierachy, x.data(), b.data(), 2, 2, cycle::wcycle, smoother::gauss_seidel, control); + // int cycles = linalg::amg_solve(hierachy, x.data(), b.data(), 2, 2, cycle::wcycle, smoother::SOR, control); // // Print solution // std::cout << "x" << std::endl; @@ -63,4 +63,4 @@ int main() // std::cout << "" << std::endl; return 0; -} \ No newline at end of file +} diff --git a/clients/examples/saamg_example.cpp b/clients/examples/saamg_example.cpp index 017d3fd..cac8ce3 100644 --- a/clients/examples/saamg_example.cpp +++ b/clients/examples/saamg_example.cpp @@ -50,9 +50,10 @@ int main() control.rel_tol = 1e-08; control.abs_tol = 1e-08; - int cycles = linalg::amg_solve(hierachy, x, b, 2, 1, linalg::Cycle::Vcycle, linalg::Smoother::Gauss_Seidel, control); - //int cycles = linalg::amg_solve(hierachy, x, b, 2, 1, linalg::Cycle::Wcycle, linalg::Smoother::Gauss_Seidel, control); - //int cycles = linalg::amg_solve(hierachy, x, b, 2, 1, linalg::Cycle::Wcycle, linalg::Smoother::Gauss_Seidel, control); + int cycles = linalg::amg_solve( + hierachy, x, b, 2, 1, linalg::cycle::vcycle, linalg::smoother::gauss_seidel, control); + //int cycles = linalg::amg_solve(hierachy, x, b, 2, 1, linalg::cycle::wcycle, linalg::smoother::gauss_seidel, control); + //int cycles = linalg::amg_solve(hierachy, x, b, 2, 1, linalg::cycle::wcycle, linalg::smoother::gauss_seidel, control); // // Print solution // std::cout << "x" << std::endl; @@ -63,4 +64,4 @@ int main() // std::cout << "" << std::endl; return 0; -} \ No newline at end of file +} diff --git a/clients/examples/uaamg_example.cpp b/clients/examples/uaamg_example.cpp index 78f0a75..4e5c1ff 100644 --- a/clients/examples/uaamg_example.cpp +++ b/clients/examples/uaamg_example.cpp @@ -50,7 +50,8 @@ int main() control.rel_tol = 1e-08; control.abs_tol = 1e-08; - int cycles = linalg::amg_solve(hierachy, x, b, 2, 2, linalg::Cycle::Wcycle, linalg::Smoother::Gauss_Seidel, control); + int cycles = linalg::amg_solve( + hierachy, x, b, 2, 2, linalg::cycle::wcycle, linalg::smoother::gauss_seidel, control); // // Print solution // std::cout << "x" << std::endl; @@ -61,4 +62,4 @@ int main() // std::cout << "" << std::endl; return 0; -} \ No newline at end of file +} diff --git a/clients/testing/CMakeLists.txt b/clients/testing/CMakeLists.txt index 9f75163..8e8632e 100644 --- a/clients/testing/CMakeLists.txt +++ b/clients/testing/CMakeLists.txt @@ -93,6 +93,7 @@ target_sources(test_main test_functions_transpose.cpp test_functions_transpose_dense.cpp test_functions_ruiz_scaling.cpp + test_functions_symmetric_ruiz_scaling.cpp test_functions_csric0.cpp test_functions_csrilu0.cpp test_functions_tridiagonal_solver.cpp diff --git a/clients/testing/test_arguments.h b/clients/testing/test_arguments.h index 5c75efb..90b7b53 100644 --- a/clients/testing/test_arguments.h +++ b/clients/testing/test_arguments.h @@ -38,11 +38,11 @@ namespace testing struct Arguments { testing::category category; // iterative_solvers, math, primitive - testing::fixture fixture; // Jacobi, CG, SpMV, exclusive_scan, etc + testing::fixture fixture; // jacobi, CG, SpMV, exclusive_scan, etc std::string group; // small, medium, large, etc std::string filename; // bmwcra_1.mtx, shipsec1.mtx, etc testing::backend backend; // CPU, GPU - testing::uplo uplo; // Lower, Upper + testing::uplo uplo; // lower, upper testing::preconditioner precond_type; testing::cycle_type cycle_type; testing::smoother_type smoother_type; @@ -62,23 +62,23 @@ namespace testing { name += "_" + backend_to_string(this->backend); } - if(this->uplo != uplo::Lower) + if(this->uplo != uplo::lower) { name += "_" + uplo_to_string(this->uplo); } - if(this->precond_type != preconditioner::None) + if(this->precond_type != preconditioner::none) { name += "_" + preconditioner_to_string(this->precond_type); } - if(this->cycle_type != cycle_type::None) + if(this->cycle_type != cycle_type::none) { name += "_" + cycle_type_to_string(this->cycle_type); } - if(this->smoother_type != smoother_type::None) + if(this->smoother_type != smoother_type::none) { name += "_" + smoother_type_to_string(this->smoother_type); } - if(this->pivoting_strategy != pivoting_strategy::None) + if(this->pivoting_strategy != pivoting_strategy::none) { name += "_" + pivoting_strategy_to_string(this->pivoting_strategy); } diff --git a/clients/testing/test_enums.h b/clients/testing/test_enums.h index 62af2ab..3f8d663 100644 --- a/clients/testing/test_enums.h +++ b/clients/testing/test_enums.h @@ -39,15 +39,15 @@ namespace testing math, primitive, csr_matrix, - Unknown + unknown }; enum class fixture { - Jacobi, - GaussSeidel, + jacobi, + gauss_seidel, SOR, - SymmGaussSeidel, + symmetric_gauss_seidel, SSOR, CG, BICGSTAB, @@ -70,7 +70,8 @@ namespace testing tridiagonal_solver, exclusive_scan, ruiz_scaling, - Unknown + symmetric_ruiz_scaling, + unknown }; enum class backend @@ -81,16 +82,16 @@ namespace testing enum class uplo { - Lower, - Upper + lower, + upper }; enum class classical_solver { - Jacobi, - GaussSeidel, + jacobi, + gauss_seidel, SOR, - SymmGaussSeidel, + symmetric_gauss_seidel, SSOR }; @@ -110,38 +111,38 @@ namespace testing enum class preconditioner { - Jacobi, - GaussSeidel, + jacobi, + gauss_seidel, SOR, - SymmGaussSeidel, + symmetric_gauss_seidel, SSOR, ILU, IC, - None + none }; enum class cycle_type { - Vcycle, - Wcycle, - Fcycle, - None + vcycle, + wcycle, + fcycle, + none }; enum class smoother_type { - Jacobi, - Gauss_Seidel, - Symm_Gauss_Seidel, + jacobi, + gauss_seidel, + symmetric_gauss_seidel, SOR, SSOR, - None + none }; enum class pivoting_strategy { - Partial, - None + partial, + none }; inline std::string category_to_string(category category) @@ -160,21 +161,21 @@ namespace testing return "csr_matrix"; } - return "Invalid"; + return "invalid"; } inline std::string fixture_to_string(fixture fixture) { switch(fixture) { - case fixture::Jacobi: - return "Jacobi"; - case fixture::GaussSeidel: - return "GaussSeidel"; + case fixture::jacobi: + return "jacobi"; + case fixture::gauss_seidel: + return "gauss_seidel"; case fixture::SOR: return "SOR"; - case fixture::SymmGaussSeidel: - return "SymmGaussSeidel"; + case fixture::symmetric_gauss_seidel: + return "symmetric_gauss_seidel"; case fixture::SSOR: return "SSOR"; case fixture::CG: @@ -219,9 +220,11 @@ namespace testing return "exclusive_scan"; case fixture::ruiz_scaling: return "ruiz_scaling"; + case fixture::symmetric_ruiz_scaling: + return "symmetric_ruiz_scaling"; } - return "Invalid"; + return "invalid"; } inline std::string backend_to_string(backend backend) @@ -234,39 +237,39 @@ namespace testing return "GPU"; } - return "Invalid"; + return "invalid"; } inline std::string uplo_to_string(uplo uplo) { switch(uplo) { - case uplo::Lower: - return "Lower"; - case uplo::Upper: - return "Upper"; + case uplo::lower: + return "lower"; + case uplo::upper: + return "upper"; } - return "Invalid"; + return "invalid"; } inline std::string classical_solver_to_string(classical_solver solver) { switch(solver) { - case classical_solver::Jacobi: - return "Jacobi"; - case classical_solver::GaussSeidel: - return "GaussSeidel"; + case classical_solver::jacobi: + return "jacobi"; + case classical_solver::gauss_seidel: + return "gauss_seidel"; case classical_solver::SOR: return "SOR"; - case classical_solver::SymmGaussSeidel: - return "SymmGaussSeidel"; + case classical_solver::symmetric_gauss_seidel: + return "symmetric_gauss_seidel"; case classical_solver::SSOR: return "SSOR"; } - return "Invalid"; + return "invalid"; } inline std::string krylov_solver_to_string(krylov_solver solver) @@ -281,7 +284,7 @@ namespace testing return "GMRES"; } - return "Invalid"; + return "invalid"; } inline std::string amg_solver_to_string(AMG_solver solver) @@ -296,83 +299,83 @@ namespace testing return "RSAMG"; } - return "Invalid"; + return "invalid"; } inline std::string preconditioner_to_string(preconditioner precond) { switch(precond) { - case preconditioner::Jacobi: - return "Jacobi"; - case preconditioner::GaussSeidel: - return "GaussSeidel"; + case preconditioner::jacobi: + return "jacobi"; + case preconditioner::gauss_seidel: + return "gauss_seidel"; case preconditioner::SOR: return "SOR"; - case preconditioner::SymmGaussSeidel: - return "SymmGaussSeidel"; + case preconditioner::symmetric_gauss_seidel: + return "symmetric_gauss_seidel"; case preconditioner::SSOR: return "SSOR"; case preconditioner::ILU: return "ILU"; case preconditioner::IC: return "IC"; - case preconditioner::None: - return "None"; + case preconditioner::none: + return "none"; } - return "Invalid"; + return "invalid"; } inline std::string cycle_type_to_string(cycle_type cycle) { switch(cycle) { - case cycle_type::Vcycle: - return "Vcycle"; - case cycle_type::Wcycle: - return "Wcycle"; - case cycle_type::Fcycle: - return "Fcycle"; - case cycle_type::None: - return "None"; + case cycle_type::vcycle: + return "vcycle"; + case cycle_type::wcycle: + return "wcycle"; + case cycle_type::fcycle: + return "fcycle"; + case cycle_type::none: + return "none"; } - return "Invalid"; + return "invalid"; } inline std::string smoother_type_to_string(smoother_type smoother) { switch(smoother) { - case smoother_type::Jacobi: - return "Jacobi"; - case smoother_type::Gauss_Seidel: - return "Gauss_Seidel"; - case smoother_type::Symm_Gauss_Seidel: - return "Symm_Gauss_Seidel"; + case smoother_type::jacobi: + return "jacobi"; + case smoother_type::gauss_seidel: + return "gauss_seidel"; + case smoother_type::symmetric_gauss_seidel: + return "symmetric_gauss_seidel"; case smoother_type::SOR: return "SOR"; case smoother_type::SSOR: return "SSOR"; - case smoother_type::None: - return "None"; + case smoother_type::none: + return "none"; } - return "Invalid"; + return "invalid"; } inline std::string pivoting_strategy_to_string(pivoting_strategy strategy) { switch(strategy) { - case pivoting_strategy::Partial: - return "Partial"; - case pivoting_strategy::None: - return "None"; + case pivoting_strategy::partial: + return "partial"; + case pivoting_strategy::none: + return "none"; } - return "Invalid"; + return "invalid"; } } diff --git a/clients/testing/test_functions.h b/clients/testing/test_functions.h index 1561332..7fe365a 100644 --- a/clients/testing/test_functions.h +++ b/clients/testing/test_functions.h @@ -52,6 +52,7 @@ namespace testing // csr matrix bool test_ruiz_scaling(Arguments arg); + bool test_symmetric_ruiz_scaling(Arguments arg); bool test_transpose(Arguments arg); bool test_transpose_dense(Arguments arg); bool test_multiply_by_vector(Arguments arg); diff --git a/clients/testing/test_functions_amg.cpp b/clients/testing/test_functions_amg.cpp index 4a6af60..68be15e 100644 --- a/clients/testing/test_functions_amg.cpp +++ b/clients/testing/test_functions_amg.cpp @@ -77,37 +77,37 @@ bool testing::test_amg(AMG_solver solver_type, Arguments arg) break; } - linalg::Cycle cycle = linalg::Cycle::Vcycle; + linalg::cycle cycle = linalg::cycle::vcycle; switch(arg.cycle_type) { - case testing::cycle_type::Vcycle: - cycle = linalg::Cycle::Vcycle; + case testing::cycle_type::vcycle: + cycle = linalg::cycle::vcycle; break; - case testing::cycle_type::Wcycle: - cycle = linalg::Cycle::Wcycle; + case testing::cycle_type::wcycle: + cycle = linalg::cycle::wcycle; break; - case testing::cycle_type::Fcycle: - cycle = linalg::Cycle::Fcycle; + case testing::cycle_type::fcycle: + cycle = linalg::cycle::fcycle; break; } - linalg::Smoother smoother = linalg::Smoother::Jacobi; + linalg::smoother smoother = linalg::smoother::jacobi; switch(arg.smoother_type) { - case testing::smoother_type::Jacobi: - smoother = linalg::Smoother::Jacobi; + case testing::smoother_type::jacobi: + smoother = linalg::smoother::jacobi; break; - case testing::smoother_type::Gauss_Seidel: - smoother = linalg::Smoother::Gauss_Seidel; + case testing::smoother_type::gauss_seidel: + smoother = linalg::smoother::gauss_seidel; break; - case testing::smoother_type::Symm_Gauss_Seidel: - smoother = linalg::Smoother::Symm_Gauss_Seidel; + case testing::smoother_type::symmetric_gauss_seidel: + smoother = linalg::smoother::symmetric_gauss_seidel; break; case testing::smoother_type::SOR: - smoother = linalg::Smoother::SOR; + smoother = linalg::smoother::SOR; break; case testing::smoother_type::SSOR: - smoother = linalg::Smoother::SSOR; + smoother = linalg::smoother::SSOR; break; } diff --git a/clients/testing/test_functions_classical.cpp b/clients/testing/test_functions_classical.cpp index 2f0a8f5..c8d5ba4 100644 --- a/clients/testing/test_functions_classical.cpp +++ b/clients/testing/test_functions_classical.cpp @@ -68,16 +68,16 @@ bool testing::test_classical(classical_solver solver_type, Arguments arg) switch(solver_type) { - case classical_solver::Jacobi: + case classical_solver::jacobi: jac.build(mat_A); break; - case classical_solver::GaussSeidel: + case classical_solver::gauss_seidel: gs.build(mat_A); break; case classical_solver::SOR: sor.build(mat_A); break; - case classical_solver::SymmGaussSeidel: + case classical_solver::symmetric_gauss_seidel: sgs.build(mat_A); break; case classical_solver::SSOR: @@ -91,16 +91,16 @@ bool testing::test_classical(classical_solver solver_type, Arguments arg) switch(solver_type) { - case classical_solver::Jacobi: + case classical_solver::jacobi: iter = jac.solve(mat_A, vec_x, vec_b, control); break; - case classical_solver::GaussSeidel: + case classical_solver::gauss_seidel: iter = gs.solve(mat_A, vec_x, vec_b, control); break; case classical_solver::SOR: iter = sor.solve(mat_A, vec_x, vec_b, control, 0.666667); break; - case classical_solver::SymmGaussSeidel: + case classical_solver::symmetric_gauss_seidel: iter = sgs.solve(mat_A, vec_x, vec_b, control); break; case classical_solver::SSOR: diff --git a/clients/testing/test_functions_dispatch.cpp b/clients/testing/test_functions_dispatch.cpp index a3b141f..8c4ef16 100644 --- a/clients/testing/test_functions_dispatch.cpp +++ b/clients/testing/test_functions_dispatch.cpp @@ -41,14 +41,14 @@ namespace testing { switch(arg.fixture) { - case fixture::Jacobi: - return test_classical(classical_solver::Jacobi, arg); - case fixture::GaussSeidel: - return test_classical(classical_solver::GaussSeidel, arg); + case fixture::jacobi: + return test_classical(classical_solver::jacobi, arg); + case fixture::gauss_seidel: + return test_classical(classical_solver::gauss_seidel, arg); case fixture::SOR: return test_classical(classical_solver::SOR, arg); - case fixture::SymmGaussSeidel: - return test_classical(classical_solver::SymmGaussSeidel, arg); + case fixture::symmetric_gauss_seidel: + return test_classical(classical_solver::symmetric_gauss_seidel, arg); case fixture::SSOR: return test_classical(classical_solver::SSOR, arg); case fixture::CG: @@ -129,6 +129,8 @@ namespace testing return test_transpose_dense(arg); case fixture::ruiz_scaling: return test_ruiz_scaling(arg); + case fixture::symmetric_ruiz_scaling: + return test_symmetric_ruiz_scaling(arg); } return false; } diff --git a/clients/testing/test_functions_krylov.cpp b/clients/testing/test_functions_krylov.cpp index 17ad9ff..27c625c 100644 --- a/clients/testing/test_functions_krylov.cpp +++ b/clients/testing/test_functions_krylov.cpp @@ -74,16 +74,16 @@ bool testing::test_krylov(krylov_solver solver_type, Arguments arg) linalg::preconditioner* p = nullptr; switch(arg.precond_type) { - case testing::preconditioner::Jacobi: + case testing::preconditioner::jacobi: p = new linalg::jacobi_precond; break; - case testing::preconditioner::GaussSeidel: + case testing::preconditioner::gauss_seidel: p = new linalg::gauss_seidel_precond; break; case testing::preconditioner::SOR: p = new linalg::SOR_precond(0.3); break; - case testing::preconditioner::SymmGaussSeidel: + case testing::preconditioner::symmetric_gauss_seidel: p = new linalg::symmetric_gauss_seidel_precond; break; case testing::preconditioner::SSOR: diff --git a/clients/testing/test_functions_ruiz_scaling.cpp b/clients/testing/test_functions_ruiz_scaling.cpp index 5cd7084..387b9f2 100644 --- a/clients/testing/test_functions_ruiz_scaling.cpp +++ b/clients/testing/test_functions_ruiz_scaling.cpp @@ -43,35 +43,18 @@ bool testing::test_ruiz_scaling(Arguments arg) linalg::csr_matrix mat_A; mat_A.read_mtx(arg.filename); - linalg::csr_matrix mat_A_copy; - mat_A_copy.copy_from(mat_A); - linalg::vector D1(mat_A.get_m()); linalg::vector D2(mat_A.get_m()); if(arg.backend == backend::GPU) { mat_A.move_to_device(); - mat_A_copy.move_to_device(); D1.move_to_device(); D2.move_to_device(); } - // Warmup - for(int i = 0; i < 4; i++) - { - mat_A.apply_ruiz_scaling(D1, D2); - mat_A.copy_from(mat_A_copy); - } - linalg::synchronize(); - - // Timed run auto t1 = std::chrono::high_resolution_clock::now(); - for(int i = 0; i < 100; i++) - { - mat_A.apply_ruiz_scaling(D1, D2); - mat_A.copy_from(mat_A_copy); - } + mat_A.apply_ruiz_scaling(D1, D2, arg.max_iters, arg.tol); linalg::synchronize(); auto t2 = std::chrono::high_resolution_clock::now(); @@ -85,5 +68,32 @@ bool testing::test_ruiz_scaling(Arguments arg) D2.move_to_host(); } + // Verify that the scaled matrix has row and column norms close to 1 + std::vector row_norm(mat_A.get_m(), 0.0); + std::vector col_norm(mat_A.get_m(), 0.0); + + for(int i = 0; i < mat_A.get_m(); i++) + { + const int start = mat_A.get_row_ptr()[i]; + const int end = mat_A.get_row_ptr()[i + 1]; + + for(int j = start; j < end; j++) + { + const int col = mat_A.get_col_ind()[j]; + const double val = mat_A.get_val()[j]; + + row_norm[i] = std::max(row_norm[i], std::abs(val)); + col_norm[col] = std::max(col_norm[col], std::abs(val)); + } + } + + for(int i = 0; i < mat_A.get_m(); i++) + { + if(std::abs(1.0 - row_norm[i]) > arg.tol || std::abs(1.0 - col_norm[i]) > arg.tol) + { + return false; + } + } + return true; } diff --git a/clients/testing/test_functions_symmetric_ruiz_scaling.cpp b/clients/testing/test_functions_symmetric_ruiz_scaling.cpp new file mode 100644 index 0000000..a268d1c --- /dev/null +++ b/clients/testing/test_functions_symmetric_ruiz_scaling.cpp @@ -0,0 +1,96 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include "test_functions.h" +#include "utility.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "linalg.h" + +bool testing::test_symmetric_ruiz_scaling(Arguments arg) +{ + std::cout << "max_iters: " << arg.max_iters << " tol: " << arg.tol << std::endl; + + linalg::csr_matrix mat_A; + mat_A.read_mtx(arg.filename); + + linalg::vector D(mat_A.get_m()); + + if(arg.backend == backend::GPU) + { + mat_A.move_to_device(); + D.move_to_device(); + } + + auto t1 = std::chrono::high_resolution_clock::now(); + mat_A.apply_symmetric_ruiz_scaling(D, arg.max_iters, arg.tol); + linalg::synchronize(); + auto t2 = std::chrono::high_resolution_clock::now(); + + std::chrono::duration ms_float = t2 - t1; + std::cout << "Solve time: " << ms_float.count() << "ms" << std::endl; + + if(arg.backend == backend::GPU) + { + mat_A.move_to_host(); + D.move_to_host(); + } + + // Verify that the scaled matrix has row and column norms close to 1 + std::vector row_norm(mat_A.get_m(), 0.0); + + for(int i = 0; i < mat_A.get_m(); i++) + { + const int start = mat_A.get_row_ptr()[i]; + const int end = mat_A.get_row_ptr()[i + 1]; + + for(int j = start; j < end; j++) + { + const int col = mat_A.get_col_ind()[j]; + const double val = mat_A.get_val()[j]; + + row_norm[i] = std::max(row_norm[i], std::abs(val)); + } + } + + for(int i = 0; i < mat_A.get_m(); i++) + { + if(std::abs(1.0 - row_norm[i]) > arg.tol) + { + return false; + } + } + + return true; +} diff --git a/clients/testing/test_functions_triangular_solve.cpp b/clients/testing/test_functions_triangular_solve.cpp index c1d30a4..253e7af 100644 --- a/clients/testing/test_functions_triangular_solve.cpp +++ b/clients/testing/test_functions_triangular_solve.cpp @@ -50,10 +50,10 @@ bool testing::test_triangular_solve(Arguments arg) switch(arg.uplo) { - case testing::uplo::Lower: + case testing::uplo::lower: mat_B.copy_lower_triangular_from(mat_A, true); break; - case testing::uplo::Upper: + case testing::uplo::upper: mat_B.copy_upper_triangular_from(mat_A, true); break; } @@ -82,10 +82,10 @@ bool testing::test_triangular_solve(Arguments arg) { switch(arg.uplo) { - case testing::uplo::Lower: + case testing::uplo::lower: mat_B.triangular_solve_lower(vec_x, vec_y, false); break; - case testing::uplo::Upper: + case testing::uplo::upper: mat_B.triangular_solve_upper(vec_x, vec_y, false); break; } @@ -98,10 +98,10 @@ bool testing::test_triangular_solve(Arguments arg) { switch(arg.uplo) { - case testing::uplo::Lower: + case testing::uplo::lower: mat_B.triangular_solve_lower(vec_x, vec_y, false); break; - case testing::uplo::Upper: + case testing::uplo::upper: mat_B.triangular_solve_upper(vec_x, vec_y, false); break; } diff --git a/clients/testing/test_functions_tridiagonal_solver.cpp b/clients/testing/test_functions_tridiagonal_solver.cpp index 4f03c6f..c517a15 100644 --- a/clients/testing/test_functions_tridiagonal_solver.cpp +++ b/clients/testing/test_functions_tridiagonal_solver.cpp @@ -88,10 +88,10 @@ bool testing::test_tridiagonal_solver(Arguments arg) linalg::pivoting_strategy pivoting; switch(arg.pivoting_strategy) { - case testing::pivoting_strategy::None: + case testing::pivoting_strategy::none: pivoting = linalg::pivoting_strategy::none; break; - case testing::pivoting_strategy::Partial: + case testing::pivoting_strategy::partial: pivoting = linalg::pivoting_strategy::partial; break; } diff --git a/clients/testing/test_yaml_loader.h b/clients/testing/test_yaml_loader.h index f6e3233..880b6e8 100644 --- a/clients/testing/test_yaml_loader.h +++ b/clients/testing/test_yaml_loader.h @@ -83,11 +83,11 @@ namespace YAML Node node; switch(rhs) { - case testing::uplo::Lower: - node = "Lower"; + case testing::uplo::lower: + node = "lower"; break; - case testing::uplo::Upper: - node = "Upper"; + case testing::uplo::upper: + node = "upper"; break; } @@ -97,13 +97,13 @@ namespace YAML static bool decode(const Node& node, testing::uplo& rhs) { std::string type = node.as(); - if(type == "Lower") + if(type == "lower") { - rhs = testing::uplo::Lower; + rhs = testing::uplo::lower; } - else if(type == "Upper") + else if(type == "upper") { - rhs = testing::uplo::Upper; + rhs = testing::uplo::upper; } return true; @@ -118,17 +118,17 @@ namespace YAML Node node; switch(rhs) { - case testing::preconditioner::Jacobi: - node = "Jacobi"; + case testing::preconditioner::jacobi: + node = "jacobi"; break; - case testing::preconditioner::GaussSeidel: - node = "GaussSeidel"; + case testing::preconditioner::gauss_seidel: + node = "gauss_seidel"; break; case testing::preconditioner::SOR: node = "SOR"; break; - case testing::preconditioner::SymmGaussSeidel: - node = "SymmGaussSeidel"; + case testing::preconditioner::symmetric_gauss_seidel: + node = "symmetric_gauss_seidel"; break; case testing::preconditioner::IC: node = "IC"; @@ -136,8 +136,8 @@ namespace YAML case testing::preconditioner::ILU: node = "ILU"; break; - case testing::preconditioner::None: - node = "None"; + case testing::preconditioner::none: + node = "none"; break; } @@ -147,21 +147,21 @@ namespace YAML static bool decode(const Node& node, testing::preconditioner& rhs) { std::string type = node.as(); - if(type == "Jacobi") + if(type == "jacobi") { - rhs = testing::preconditioner::Jacobi; + rhs = testing::preconditioner::jacobi; } - if(type == "GaussSeidel") + if(type == "gauss_seidel") { - rhs = testing::preconditioner::GaussSeidel; + rhs = testing::preconditioner::gauss_seidel; } if(type == "SOR") { rhs = testing::preconditioner::SOR; } - if(type == "SymmGaussSeidel") + if(type == "symmetric_gauss_seidel") { - rhs = testing::preconditioner::SymmGaussSeidel; + rhs = testing::preconditioner::symmetric_gauss_seidel; } if(type == "SSOR") { @@ -175,9 +175,9 @@ namespace YAML { rhs = testing::preconditioner::ILU; } - else if(type == "None") + else if(type == "none") { - rhs = testing::preconditioner::None; + rhs = testing::preconditioner::none; } return true; @@ -192,14 +192,14 @@ namespace YAML Node node; switch(rhs) { - case testing::cycle_type::Fcycle: - node = "Fcycle"; + case testing::cycle_type::fcycle: + node = "fcycle"; break; - case testing::cycle_type::Vcycle: - node = "Vcycle"; + case testing::cycle_type::vcycle: + node = "vcycle"; break; - case testing::cycle_type::Wcycle: - node = "Wcycle"; + case testing::cycle_type::wcycle: + node = "wcycle"; break; } @@ -209,17 +209,17 @@ namespace YAML static bool decode(const Node& node, testing::cycle_type& rhs) { std::string type = node.as(); - if(type == "Fcycle") + if(type == "fcycle") { - rhs = testing::cycle_type::Fcycle; + rhs = testing::cycle_type::fcycle; } - else if(type == "Vcycle") + else if(type == "vcycle") { - rhs = testing::cycle_type::Vcycle; + rhs = testing::cycle_type::vcycle; } - else if(type == "Wcycle") + else if(type == "wcycle") { - rhs = testing::cycle_type::Wcycle; + rhs = testing::cycle_type::wcycle; } return true; @@ -234,14 +234,14 @@ namespace YAML Node node; switch(rhs) { - case testing::smoother_type::Jacobi: - node = "Jacobi"; + case testing::smoother_type::jacobi: + node = "jacobi"; break; - case testing::smoother_type::Gauss_Seidel: - node = "Gauss_Seidel"; + case testing::smoother_type::gauss_seidel: + node = "gauss_seidel"; break; - case testing::smoother_type::Symm_Gauss_Seidel: - node = "Symm_Gauss_Seidel"; + case testing::smoother_type::symmetric_gauss_seidel: + node = "symmetric_gauss_seidel"; break; case testing::smoother_type::SOR: node = "SOR"; @@ -257,17 +257,17 @@ namespace YAML static bool decode(const Node& node, testing::smoother_type& rhs) { std::string type = node.as(); - if(type == "Jacobi") + if(type == "jacobi") { - rhs = testing::smoother_type::Jacobi; + rhs = testing::smoother_type::jacobi; } - else if(type == "Gauss_Seidel") + else if(type == "gauss_seidel") { - rhs = testing::smoother_type::Gauss_Seidel; + rhs = testing::smoother_type::gauss_seidel; } - else if(type == "Symm_Gauss_Seidel") + else if(type == "symmetric_gauss_seidel") { - rhs = testing::smoother_type::Symm_Gauss_Seidel; + rhs = testing::smoother_type::symmetric_gauss_seidel; } else if(type == "SOR") { @@ -290,11 +290,11 @@ namespace YAML Node node; switch(rhs) { - case testing::pivoting_strategy::Partial: - node = "Partial"; + case testing::pivoting_strategy::partial: + node = "partial"; break; - case testing::pivoting_strategy::None: - node = "None"; + case testing::pivoting_strategy::none: + node = "none"; break; } @@ -304,13 +304,13 @@ namespace YAML static bool decode(const Node& node, testing::pivoting_strategy& rhs) { std::string type = node.as(); - if(type == "Partial") + if(type == "partial") { - rhs = testing::pivoting_strategy::Partial; + rhs = testing::pivoting_strategy::partial; } - else if(type == "None") + else if(type == "none") { - rhs = testing::pivoting_strategy::None; + rhs = testing::pivoting_strategy::none; } return true; @@ -373,17 +373,17 @@ inline testing::category string_to_category(const std::string& str) return it->second; } - return testing::category::Unknown; + return testing::category::unknown; } inline testing::fixture string_to_fixture(const std::string& str) { // Static map for efficiency. It's initialized only once. static const std::unordered_map fixture_map - = {{"Jacobi", testing::fixture::Jacobi}, - {"GaussSeidel", testing::fixture::GaussSeidel}, + = {{"jacobi", testing::fixture::jacobi}, + {"gauss_seidel", testing::fixture::gauss_seidel}, {"SOR", testing::fixture::SOR}, - {"SymmGaussSeidel", testing::fixture::SymmGaussSeidel}, + {"symmetric_gauss_seidel", testing::fixture::symmetric_gauss_seidel}, {"SSOR", testing::fixture::SSOR}, {"CG", testing::fixture::CG}, {"BICGSTAB", testing::fixture::BICGSTAB}, @@ -405,6 +405,8 @@ inline testing::fixture string_to_fixture(const std::string& str) {"CSRIC0", testing::fixture::CSRIC0}, {"CSRILU0", testing::fixture::CSRILU0}, {"tridiagonal_solver", testing::fixture::tridiagonal_solver}, + {"ruiz_scaling", testing::fixture::ruiz_scaling}, + {"symmetric_ruiz_scaling", testing::fixture::symmetric_ruiz_scaling}, {"exclusive_scan", testing::fixture::exclusive_scan}}; // Find the string in the map @@ -416,7 +418,7 @@ inline testing::fixture string_to_fixture(const std::string& str) return it->second; } - return testing::fixture::Unknown; + return testing::fixture::unknown; } // Helper struct to hold all parameter vectors @@ -497,12 +499,12 @@ inline std::vector generate_tests(const std::string category TestParameters params; params.matrices = read_group_values("matrix_file", std::string("")); params.backends = read_group_values("backend", testing::backend::CPU); - params.uplos = read_group_values("uplo", testing::uplo::Lower); - params.precond_types = read_group_values("precond", testing::preconditioner::None); - params.cycle_types = read_group_values("cycle", testing::cycle_type::None); - params.smoother_types = read_group_values("smoother", testing::smoother_type::None); + params.uplos = read_group_values("uplo", testing::uplo::lower); + params.precond_types = read_group_values("precond", testing::preconditioner::none); + params.cycle_types = read_group_values("cycle", testing::cycle_type::none); + params.smoother_types = read_group_values("smoother", testing::smoother_type::none); params.pivoting_strategies - = read_group_values("pivoting_strategy", testing::pivoting_strategy::None); + = read_group_values("pivoting_strategy", testing::pivoting_strategy::none); params.presmoothings = read_group_values("presmoothing", -1); params.postsmoothings = read_group_values("postsmoothing", -1); params.max_iters = read_group_values("max_iters", -1); diff --git a/clients/testing/tests/CMakeLists.txt b/clients/testing/tests/CMakeLists.txt index a918719..58cd732 100644 --- a/clients/testing/tests/CMakeLists.txt +++ b/clients/testing/tests/CMakeLists.txt @@ -50,4 +50,5 @@ target_sources(test_main test_tridiagonal_solver.cpp test_exclusive_scan.cpp test_ruiz_scaling.cpp + test_symmetric_ruiz_scaling.cpp ) diff --git a/clients/testing/tests/test_BICGSTAB.yaml b/clients/testing/tests/test_BICGSTAB.yaml index 84da881..fa0ddcd 100644 --- a/clients/testing/tests/test_BICGSTAB.yaml +++ b/clients/testing/tests/test_BICGSTAB.yaml @@ -1,6 +1,6 @@ Tests: small: - precond: [None, Jacobi] + precond: [none, jacobi] matrix_file: ["matrices/SPD/mesh1em6/mesh1em6.mtx", "matrices/SPD/bcsstk01/bcsstk01.mtx", "matrices/SPD/mesh3em5/mesh3em5.mtx", diff --git a/clients/testing/tests/test_CG.yaml b/clients/testing/tests/test_CG.yaml index 0f3e2d8..4ca0e40 100644 --- a/clients/testing/tests/test_CG.yaml +++ b/clients/testing/tests/test_CG.yaml @@ -1,6 +1,6 @@ Tests: small1: - precond: [None, Jacobi, SOR] + precond: [none, jacobi, SOR] matrix_file: ["matrices/SPD/bcsstm02/bcsstm02.mtx", "matrices/SPD/mesh1em6/mesh1em6.mtx", "matrices/SPD/mesh3em5/mesh3em5.mtx", @@ -17,7 +17,7 @@ Tests: backend: [CPU, GPU] small2: - precond: [Jacobi] + precond: [jacobi] matrix_file: ["matrices/SPD/nos6/nos6.mtx", "matrices/SPD/bodyy6/bodyy6.mtx", "matrices/SPD/cvxbqp1/cvxbqp1.mtx", @@ -29,7 +29,7 @@ Tests: backend: [CPU, GPU] small3: - precond: [SymmGaussSeidel, SSOR] + precond: [symmetric_gauss_seidel, SSOR] matrix_file: ["matrices/SPD/bcsstm02/bcsstm02.mtx", "matrices/SPD/mesh1em6/mesh1em6.mtx", "matrices/SPD/mesh3em5/mesh3em5.mtx", @@ -47,7 +47,7 @@ Tests: backend: [CPU, GPU] medium1: - precond: [Jacobi] + precond: [jacobi] matrix_file: ["matrices/SPD/Dubcova2/Dubcova2.mtx", "matrices/SPD/denormal/denormal.mtx", "matrices/SPD/raefsky4/raefsky4.mtx", @@ -61,7 +61,7 @@ Tests: backend: [CPU, GPU] medium2: - precond: [Jacobi] + precond: [jacobi] matrix_file: ["matrices/SPD/2cubes_sphere/2cubes_sphere.mtx", "matrices/SPD/parabolic_fem/parabolic_fem.mtx", "matrices/SPD/offshore/offshore.mtx"] @@ -69,7 +69,7 @@ Tests: backend: [CPU, GPU] medium3: - precond: [SymmGaussSeidel, SSOR] + precond: [symmetric_gauss_seidel, SSOR] matrix_file: ["matrices/SPD/denormal/denormal.mtx", "matrices/SPD/raefsky4/raefsky4.mtx", "matrices/SPD/thermomech_dM/thermomech_dM.mtx", @@ -80,14 +80,14 @@ Tests: backend: [CPU, GPU] large1: - precond: [Jacobi] + precond: [jacobi] matrix_file: ["matrices/SPD/crankseg_1/crankseg_1.mtx", "matrices/SPD/crankseg_2/crankseg_2.mtx"] max_iters: [300] backend: [CPU, GPU] slow2: - precond: [Jacobi] + precond: [jacobi] matrix_file: ["matrices/SPD/nos1/nos1.mtx", "matrices/SPD/Kuu/Kuu.mtx", "matrices/SPD/s1rmq4m1/s1rmq4m1.mtx", diff --git a/clients/testing/tests/test_GMRES.yaml b/clients/testing/tests/test_GMRES.yaml index 9c5b127..64f01aa 100644 --- a/clients/testing/tests/test_GMRES.yaml +++ b/clients/testing/tests/test_GMRES.yaml @@ -1,6 +1,6 @@ Tests: small: - precond: [None, Jacobi] + precond: [none, jacobi] matrix_file: ["matrices/SPD/crystm02/crystm02.mtx", "matrices/SPD/thermomech_TK/thermomech_TK.mtx"] max_iters: [1000] diff --git a/clients/testing/tests/test_gauss_seidel.cpp b/clients/testing/tests/test_gauss_seidel.cpp index 87f0c03..feb5e69 100644 --- a/clients/testing/tests/test_gauss_seidel.cpp +++ b/clients/testing/tests/test_gauss_seidel.cpp @@ -26,4 +26,4 @@ #include "../test.h" -INSTANTIATE_TEST(iterative_solvers, GaussSeidel, classical, "tests/test_gauss_seidel.yaml"); +INSTANTIATE_TEST(iterative_solvers, gauss_seidel, classical, "tests/test_gauss_seidel.yaml"); diff --git a/clients/testing/tests/test_jacobi.cpp b/clients/testing/tests/test_jacobi.cpp index 9d668d5..7f0f3ff 100644 --- a/clients/testing/tests/test_jacobi.cpp +++ b/clients/testing/tests/test_jacobi.cpp @@ -26,4 +26,4 @@ #include "../test.h" -INSTANTIATE_TEST(iterative_solvers, Jacobi, classical, "tests/test_jacobi.yaml"); +INSTANTIATE_TEST(iterative_solvers, jacobi, classical, "tests/test_jacobi.yaml"); diff --git a/clients/testing/tests/test_ruiz_scaling.yaml b/clients/testing/tests/test_ruiz_scaling.yaml index c36d9d0..28461f9 100644 --- a/clients/testing/tests/test_ruiz_scaling.yaml +++ b/clients/testing/tests/test_ruiz_scaling.yaml @@ -26,4 +26,53 @@ Tests: "matrices/SPD/wathen100/wathen100.mtx", "matrices/SPD/apache1/apache1.mtx", "matrices/SPD/wathen120/wathen120.mtx"] + max_iters: [30] + tol: [1e-6] + backend: [CPU] + + medium: + matrix_file: ["matrices/SPD/olafu/olafu.mtx", + "matrices/SPD/gyro/gyro.mtx", + "matrices/SPD/Dubcova2/Dubcova2.mtx", + "matrices/SPD/msc23052/msc23052.mtx", + "matrices/SPD/bcsstk36/bcsstk36.mtx", + "matrices/SPD/msc10848/msc10848.mtx", + "matrices/SPD/raefsky4/raefsky4.mtx", + "matrices/SPD/thermomech_dM/thermomech_dM.mtx", + "matrices/SPD/2cubes_sphere/2cubes_sphere.mtx", + "matrices/SPD/qa8fm/qa8fm.mtx", + "matrices/SPD/denormal/denormal.mtx", + "matrices/SPD/cfd1/cfd1.mtx", + "matrices/SPD/oilpan/oilpan.mtx", + "matrices/SPD/vanbody/vanbody.mtx", + "matrices/SPD/ct20stif/ct20stif.mtx", + "matrices/SPD/cfd2/cfd2.mtx", + "matrices/SPD/shipsec8/shipsec8.mtx", + "matrices/SPD/shipsec1/shipsec1.mtx", + "matrices/SPD/Dubcova3/Dubcova3.mtx", + "matrices/SPD/parabolic_fem/parabolic_fem.mtx", + "matrices/SPD/s3dkt3m2/s3dkt3m2.mtx", + "matrices/SPD/smt/smt.mtx", + "matrices/SPD/ship_003/ship_003.mtx", + "matrices/SPD/ship_001/ship_001.mtx", + "matrices/SPD/cant/cant.mtx", + "matrices/SPD/offshore/offshore.mtx", + "matrices/SPD/pdb1HYS/pdb1HYS.mtx", + "matrices/SPD/s3dkq4m2/s3dkq4m2.mtx", + "matrices/SPD/thread/thread.mtx", + "matrices/SPD/shipsec5/shipsec5.mtx", + "matrices/SPD/apache2/apache2.mtx", + "matrices/SPD/ecology2/ecology2.mtx", + "matrices/SPD/tmt_sym/tmt_sym.mtx", + "matrices/SPD/boneS01/boneS01.mtx", + "matrices/SPD/consph/consph.mtx", + "matrices/SPD/nd6k/nd6k.mtx", + "matrices/SPD/bmw7st_1/bmw7st_1.mtx", + "matrices/SPD/G3_circuit/G3_circuit.mtx", + "matrices/SPD/thermal2/thermal2.mtx", + "matrices/SPD/x104/x104.mtx", + "matrices/SPD/m_t1/m_t1.mtx", + "matrices/SPD/hood/hood.mtx"] + max_iters: [30] + tol: [1e-6] backend: [CPU] diff --git a/clients/testing/tests/test_symmetric_gauss_seidel.cpp b/clients/testing/tests/test_symmetric_gauss_seidel.cpp index 2fd43aa..a2c191d 100644 --- a/clients/testing/tests/test_symmetric_gauss_seidel.cpp +++ b/clients/testing/tests/test_symmetric_gauss_seidel.cpp @@ -27,6 +27,6 @@ #include "../test.h" INSTANTIATE_TEST(iterative_solvers, - SymmGaussSeidel, + symmetric_gauss_seidel, classical, "tests/test_symmetric_gauss_seidel.yaml"); diff --git a/clients/testing/tests/test_symmetric_ruiz_scaling.cpp b/clients/testing/tests/test_symmetric_ruiz_scaling.cpp new file mode 100644 index 0000000..ee46cd1 --- /dev/null +++ b/clients/testing/tests/test_symmetric_ruiz_scaling.cpp @@ -0,0 +1,32 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#include "../test.h" + +INSTANTIATE_TEST(csr_matrix, + symmetric_ruiz_scaling, + symmetric_ruiz_scaling, + "tests/test_symmetric_ruiz_scaling.yaml"); diff --git a/clients/testing/tests/test_symmetric_ruiz_scaling.yaml b/clients/testing/tests/test_symmetric_ruiz_scaling.yaml new file mode 100644 index 0000000..216c43b --- /dev/null +++ b/clients/testing/tests/test_symmetric_ruiz_scaling.yaml @@ -0,0 +1,78 @@ +Tests: + small: + matrix_file: ["matrices/SPD/bcsstm02/bcsstm02.mtx", + "matrices/SPD/ex5/ex5.mtx", + "matrices/SPD/mesh1em6/mesh1em6.mtx", + "matrices/SPD/bcsstk01/bcsstk01.mtx", + "matrices/SPD/mesh3em5/mesh3em5.mtx", + "matrices/SPD/mesh2em5/mesh2em5.mtx", + "matrices/SPD/fv1/fv1.mtx", + "matrices/SPD/fv2/fv2.mtx", + "matrices/SPD/fv3/fv3.mtx", + "matrices/SPD/bodyy4/bodyy4.mtx", + "matrices/SPD/shallow_water1/shallow_water1.mtx", + "matrices/SPD/shallow_water2/shallow_water2.mtx", + "matrices/SPD/gyro_m/gyro_m.mtx", + "matrices/SPD/crystm02/crystm02.mtx", + "matrices/SPD/crystm03/crystm03.mtx", + "matrices/SPD/finan512/finan512.mtx", + "matrices/SPD/thermomech_TC/thermomech_TC.mtx", + "matrices/SPD/thermomech_TK/thermomech_TK.mtx", + "matrices/SPD/bundle1/bundle1.mtx", + "matrices/SPD/nos6/nos6.mtx", + "matrices/SPD/bodyy6/bodyy6.mtx", + "matrices/SPD/cvxbqp1/cvxbqp1.mtx", + "matrices/SPD/bcsstk38/bcsstk38.mtx", + "matrices/SPD/wathen100/wathen100.mtx", + "matrices/SPD/apache1/apache1.mtx", + "matrices/SPD/wathen120/wathen120.mtx"] + max_iters: [30] + tol: [1e-6] + backend: [GPU] + + medium: + matrix_file: ["matrices/SPD/olafu/olafu.mtx", + "matrices/SPD/gyro/gyro.mtx", + "matrices/SPD/Dubcova2/Dubcova2.mtx", + "matrices/SPD/msc23052/msc23052.mtx", + "matrices/SPD/bcsstk36/bcsstk36.mtx", + "matrices/SPD/msc10848/msc10848.mtx", + "matrices/SPD/raefsky4/raefsky4.mtx", + "matrices/SPD/thermomech_dM/thermomech_dM.mtx", + "matrices/SPD/2cubes_sphere/2cubes_sphere.mtx", + "matrices/SPD/qa8fm/qa8fm.mtx", + "matrices/SPD/denormal/denormal.mtx", + "matrices/SPD/cfd1/cfd1.mtx", + "matrices/SPD/oilpan/oilpan.mtx", + "matrices/SPD/vanbody/vanbody.mtx", + "matrices/SPD/ct20stif/ct20stif.mtx", + "matrices/SPD/cfd2/cfd2.mtx", + "matrices/SPD/shipsec8/shipsec8.mtx", + "matrices/SPD/shipsec1/shipsec1.mtx", + "matrices/SPD/Dubcova3/Dubcova3.mtx", + "matrices/SPD/parabolic_fem/parabolic_fem.mtx", + "matrices/SPD/s3dkt3m2/s3dkt3m2.mtx", + "matrices/SPD/smt/smt.mtx", + "matrices/SPD/ship_003/ship_003.mtx", + "matrices/SPD/ship_001/ship_001.mtx", + "matrices/SPD/cant/cant.mtx", + "matrices/SPD/offshore/offshore.mtx", + "matrices/SPD/pdb1HYS/pdb1HYS.mtx", + "matrices/SPD/s3dkq4m2/s3dkq4m2.mtx", + "matrices/SPD/thread/thread.mtx", + "matrices/SPD/shipsec5/shipsec5.mtx", + "matrices/SPD/apache2/apache2.mtx", + "matrices/SPD/ecology2/ecology2.mtx", + "matrices/SPD/tmt_sym/tmt_sym.mtx", + "matrices/SPD/boneS01/boneS01.mtx", + "matrices/SPD/consph/consph.mtx", + "matrices/SPD/nd6k/nd6k.mtx", + "matrices/SPD/bmw7st_1/bmw7st_1.mtx", + "matrices/SPD/G3_circuit/G3_circuit.mtx", + "matrices/SPD/thermal2/thermal2.mtx", + "matrices/SPD/x104/x104.mtx", + "matrices/SPD/m_t1/m_t1.mtx", + "matrices/SPD/hood/hood.mtx"] + max_iters: [30] + tol: [1e-6] + backend: [GPU] diff --git a/clients/testing/tests/test_triangular_solve.yaml b/clients/testing/tests/test_triangular_solve.yaml index c09c0ea..615964b 100644 --- a/clients/testing/tests/test_triangular_solve.yaml +++ b/clients/testing/tests/test_triangular_solve.yaml @@ -16,7 +16,7 @@ Tests: "matrices/SPD/thermomech_TC/thermomech_TC.mtx", "matrices/SPD/thermomech_TK/thermomech_TK.mtx", "matrices/SPD/bundle1/bundle1.mtx"] - uplo: [Lower, Upper] + uplo: [lower, upper] backend: [CPU, GPU] medium: @@ -59,5 +59,5 @@ Tests: "matrices/SPD/thermal2/thermal2.mtx", "matrices/SPD/m_t1/m_t1.mtx", "matrices/SPD/hood/hood.mtx"] - uplo: [Lower, Upper] + uplo: [lower, upper] backend: [CPU, GPU] diff --git a/clients/testing/tests/test_tridiagonal_solver.yaml b/clients/testing/tests/test_tridiagonal_solver.yaml index cbb72a8..f178392 100644 --- a/clients/testing/tests/test_tridiagonal_solver.yaml +++ b/clients/testing/tests/test_tridiagonal_solver.yaml @@ -2,41 +2,41 @@ Tests: small: m: [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] n: [2, 3, 4, 5, 6, 7, 8, 9, 10] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [CPU, GPU] small1: m: [20, 40, 75, 110, 150, 312, 625, 1000] n: [11, 12, 13, 14, 15, 16, 17, 18, 19, 20] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [CPU, GPU] small2: m: [1280, 2560, 4096, 11058, 27384, 41056, 77028, 131072] n: [21, 22, 23, 24, 25, 26, 27, 28, 29, 30] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [GPU] medium: m: [8, 29, 55, 96, 131, 175, 256, 888] n: [456, 789, 1024] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [GPU] medium1: m: [1768, 3000, 15673, 22096, 34386, 61093, 100972] n: [1200, 1350, 1587, 2000] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [GPU] large: m: [256, 512, 1024, 2048, 4096] n: [4096, 8192, 16384, 32768, 65536, 131072] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [GPU] large1: m: [130000, 155789, 475890, 1500600, 4962341, 16500700] n: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - pivoting_strategy: [None, Partial] + pivoting_strategy: [none, partial] backend: [GPU] diff --git a/library/include/csr_matrix.h b/library/include/csr_matrix.h index e8eda1a..ae299a9 100644 --- a/library/include/csr_matrix.h +++ b/library/include/csr_matrix.h @@ -365,7 +365,8 @@ namespace linalg */ void transpose(csr_matrix& T) const; - void apply_ruiz_scaling(vector& D1, vector& D2); + void apply_ruiz_scaling(vector& D1, vector& D2, int max_k, double tol); + void apply_symmetric_ruiz_scaling(vector& D, int max_k, double tol); /*! \brief Reads a sparse matrix from a Matrix Market (.mtx) file into this CSR object. * \param filename The path to the .mtx file. diff --git a/library/include/iterative_solvers/amg/amg.h b/library/include/iterative_solvers/amg/amg.h index 532c994..92c3159 100644 --- a/library/include/iterative_solvers/amg/amg.h +++ b/library/include/iterative_solvers/amg/amg.h @@ -152,7 +152,7 @@ namespace linalg * at various levels of the multigrid hierarchy. This enumeration defines the * commonly used cycle types. */ - enum class Cycle + enum class cycle { /*! \brief V-cycle. * @@ -162,7 +162,7 @@ namespace linalg * back to the finer level followed by another smoothing step. The process resembles * the shape of the letter "V". */ - Vcycle, + vcycle, /*! \brief W-cycle. * @@ -172,7 +172,7 @@ namespace linalg * can be more effective at reducing low-frequency errors but is also more computationally * expensive. The process resembles the shape of the letter "W". */ - Wcycle, + wcycle, /*! \brief F-cycle (or Full Multigrid cycle). * @@ -183,7 +183,7 @@ namespace linalg * and often leads to faster convergence overall. The process resembles the shape * of the letter "F" if visualized across multiple levels. */ - Fcycle + fcycle }; /*! \ingroup iterative_solvers_amg @@ -196,7 +196,7 @@ namespace linalg * smoother can significantly impact the convergence rate and efficiency of the * amg method. */ - enum class Smoother + enum class smoother { /*! \brief Jacobi smoother. * @@ -204,7 +204,7 @@ namespace linalg * based on the values of all other unknowns from the previous iteration. It is * easy to parallelize but often converges slowly. */ - Jacobi, + jacobi, /*! \brief Gauss-Seidel smoother. * @@ -213,7 +213,7 @@ namespace linalg * iteration. This sequential dependency can lead to faster convergence than Jacobi * in many cases. */ - Gauss_Seidel, + gauss_seidel, /*! \brief Symmetric Gauss-Seidel (SGS) smoother. * @@ -222,7 +222,7 @@ namespace linalg * application can improve the smoothing properties, particularly for symmetric * positive definite systems. */ - Symm_Gauss_Seidel, + symmetric_gauss_seidel, /*! \brief Successive Over-Relaxation (SOR) smoother. * @@ -277,12 +277,12 @@ namespace linalg * high-frequency errors. * * The `Cycle` enumeration defines how the coarse-grid problem (step 2c) is handled recursively: - * - \b V-Cycle (`Cycle::Vcycle`): The coarse-grid problem is solved by recursively applying + * - \b V-Cycle (`cycle::vcycle`): The coarse-grid problem is solved by recursively applying * a single amg V-cycle. This is the simplest and most common cycle type. - * - \b W-Cycle (`Cycle::Wcycle`): The coarse-grid problem is solved by recursively applying + * - \b W-Cycle (`cycle::wcycle`): The coarse-grid problem is solved by recursively applying * *more than one* V-cycle (typically two) on the coarser level before interpolating back. * This can lead to more thorough error reduction at the cost of increased computation per cycle. - * - \b F-Cycle (`Cycle::Fcycle`): The F-cycle starts by recursively solving the problem on the + * - \b F-Cycle (`cycle::fcycle`): The F-cycle starts by recursively solving the problem on the * coarsest level, then interpolates the solution to the next finer level and performs * one or more V-cycles. This process is repeated, gradually moving to finer levels. * It's often used to provide a good initial guess and can lead to faster overall convergence. @@ -296,9 +296,9 @@ namespace linalg * \param b The right-hand side vector of the linear system. * \param n1 The number of pre-smoothing iterations to perform at each level. * \param n2 The number of post-smoothing iterations to perform at each level. - * \param cycle An enumeration value (`Cycle::Vcycle`, `Cycle::Wcycle`, or `Cycle::Fcycle`) + * \param cycle An enumeration value (`cycle::vcycle`, `cycle::wcycle`, or `cycle::fcycle`) * specifying the type of multigrid cycle to be performed. - * \param smoother An enumeration value (`Smoother::Jacobi`, `Smoother::Gauss_Seidel`, `Smoother::Symm_Gauss_Seidel`, `Smoother::SOR`, `Smoother::SSOR`) + * \param smoother An enumeration value (`smoother::jacobi`, `smoother::gauss_seidel`, `smoother::symmetric_gauss_seidel`, `smoother::SOR`, `smoother::SSOR`) * specifying the type of stationary iterative method to be used as a smoother at each level. * \param control An `iter_control` object that manages the overall iterative process, * including convergence tolerance (`rel_tol`, `abs_tol`) and maximum number of cycles (`max_cycle`). @@ -310,8 +310,8 @@ namespace linalg * * \section amg_example Example Usage * Below is a simplified example demonstrating how to use the `amg_solve` function. - * This assumes `csr_matrix`, `vector`, `hierarchy`, `iter_control`, `Cycle` enum, - * and `Smoother` enum are properly defined and functional. + * This assumes `csr_matrix`, `vector`, `hierarchy`, `iter_control`, `cycle` enum, + * and `smoother` enum are properly defined and functional. * * \code * #include "linalglib.h" @@ -379,8 +379,8 @@ namespace linalg * // 3. Set amg parameters * int n1 = 2; // Number of pre-smoothing iterations * int n2 = 2; // Number of post-smoothing iterations - * Cycle cycle_type = Cycle::Vcycle; // Use V-cycle - * Smoother smoother_type = Smoother::Jacobi; // Use Jacobi smoother + * cycle cycle_type = cycle::vcycle; // Use V-cycle + * smoother smoother_type = smoother::jacobi; // Use Jacobi smoother * * // 4. Set up iteration control * iter_control control; @@ -413,8 +413,8 @@ namespace linalg const vector& b, int n1, int n2, - Cycle cycle, - Smoother smoother, + cycle cycle, + smoother smoother, iter_control control); } diff --git a/library/src/backend/device/cuda/cuda_ruiz_scaling.cu b/library/src/backend/device/cuda/cuda_ruiz_scaling.cu index 06adb45..86aee49 100644 --- a/library/src/backend/device/cuda/cuda_ruiz_scaling.cu +++ b/library/src/backend/device/cuda/cuda_ruiz_scaling.cu @@ -28,19 +28,99 @@ #include "cuda_ruiz_scaling.h" +#include "ruiz_scaling_kernels.cuh" + #include "../../../trace.h" //------------------------------------------------------------------------------- // Ruiz scaling //------------------------------------------------------------------------------- template -void linalg::cuda_ruiz_scaling(T* D1, - const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, T* D2) +void linalg::cuda_ruiz_scaling(T* D1, + const int* csr_row_ptr, + const int* csr_col_ind, + T* csr_val, + int m, + T* D2, + int max_k, + T tol) { ROUTINE_TRACE("linalg::cuda_ruiz_scaling_impl"); - //scale_diagonal_kernel<256> - // <<<((m - 1) / 256 + 1), 256>>>(m, csr_row_ptr, csr_col_ind, csr_val, scalar); - //CHECK_CUDA_LAUNCH_ERROR(); } -template void linalg::cuda_ruiz_scaling(double* D1, const int* csr_row_ptr, const int* csr_col_ind, double* csr_val, int m, double* D2); +template +void linalg::cuda_symmetric_ruiz_scaling( + T* D, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, int max_k, T tol) +{ + ROUTINE_TRACE("linalg::cuda_symmetric_ruiz_scaling_impl"); + + T* workspace = nullptr; + CHECK_CUDA(cudaMalloc((void**)&workspace, sizeof(T) * 256)); + CHECK_CUDA(cudaMemset((void**)&workspace, 0, sizeof(T) * 256)); + + T* dmax_divergence = nullptr; + CHECK_CUDA(cudaMalloc((void**)&dmax_divergence, sizeof(T))); + + // D^0 = I + set_D_to_one_kernel<256><<<((m - 1) / 256 + 1), 256>>>(m, D); + CHECK_CUDA_LAUNCH_ERROR(); + + // DR = DC + T* DR = nullptr; + CHECK_CUDA(cudaMalloc((void**)&DR, sizeof(T) * m)); + + for(int iter = 0; iter < max_k; iter++) + { + CHECK_CUDA(cudaMemset((void**)&DR, 0, sizeof(T) * m)); + + fill_DR_kernel<256, 32><<<((m - 1) / (256 / 32) + 1), 256>>>(m, csr_row_ptr, csr_val, DR); + CHECK_CUDA_LAUNCH_ERROR(); + + const T eps = std::numeric_limits::epsilon(); + + compute_max_divergence_part1<256><<<256, 256>>>(m, eps, DR, workspace); + CHECK_CUDA_LAUNCH_ERROR(); + + compute_max_divergence_part2<256><<<1, 256>>>(m, workspace, dmax_divergence); + CHECK_CUDA_LAUNCH_ERROR(); + + T max_divergence = static_cast(0); + CHECK_CUDA(cudaMemcpy(&max_divergence, dmax_divergence, sizeof(T), cudaMemcpyDeviceToHost)); + + std::cout << "max_divergence: " << max_divergence << std::endl; + + if(max_divergence < tol) + { + break; + } + + // Update A_k+1 = DR^-1 * A_k * DR^-1 + update_A<256, 32> + <<<((m - 1) / (256 / 32) + 1), 256>>>(m, csr_row_ptr, csr_col_ind, csr_val, DR); + CHECK_CUDA_LAUNCH_ERROR(); + + // Update D_k+1 = D_k * DR^-1 + update_D<256><<<((m - 1) / 256 + 1), 256>>>(m, DR, D); + CHECK_CUDA_LAUNCH_ERROR(); + } + + CHECK_CUDA(cudaFree(DR)); + CHECK_CUDA(cudaFree(workspace)); + CHECK_CUDA(cudaFree(dmax_divergence)); +} + +template void linalg::cuda_ruiz_scaling(double* D1, + const int* csr_row_ptr, + const int* csr_col_ind, + double* csr_val, + int m, + double* D2, + int max_k, + double tol); +template void linalg::cuda_symmetric_ruiz_scaling(double* D, + const int* csr_row_ptr, + const int* csr_col_ind, + double* csr_val, + int m, + int max_k, + double tol); diff --git a/library/src/backend/device/cuda/cuda_ruiz_scaling.h b/library/src/backend/device/cuda/cuda_ruiz_scaling.h index 8a447ae..8dd9d46 100644 --- a/library/src/backend/device/cuda/cuda_ruiz_scaling.h +++ b/library/src/backend/device/cuda/cuda_ruiz_scaling.h @@ -29,8 +29,18 @@ namespace linalg { template - void cuda_ruiz_scaling( - T* D1, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, T* D2); + void cuda_ruiz_scaling(T* D1, + const int* csr_row_ptr, + const int* csr_col_ind, + T* csr_val, + int m, + T* D2, + int max_k, + T tol); + + template + void cuda_symmetric_ruiz_scaling( + T* D, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, int max_k, T tol); } #endif diff --git a/library/src/backend/device/cuda/ruiz_scaling_kernels.cuh b/library/src/backend/device/cuda/ruiz_scaling_kernels.cuh new file mode 100644 index 0000000..41197ef --- /dev/null +++ b/library/src/backend/device/cuda/ruiz_scaling_kernels.cuh @@ -0,0 +1,185 @@ +//******************************************************************************** +// +// MIT License +// +// Copyright(c) 2026 James Sandham +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this softwareand associated documentation files(the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and /or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions : +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +//******************************************************************************** + +#ifndef RUIZ_SCALING_KERNELS_KERNELS_H +#define RUIZ_SCALING_KERNELS_KERNELS_H + +#include "common.cuh" + +template +__global__ void set_D_to_one_kernel(int m, T* __restrict__ DR) +{ + int tid = threadIdx.x; + int bid = blockIdx.x; + int gid = tid + BLOCKSIZE * bid; + + if(gid < m) + { + DR[gid] = static_cast(1); + } +} + +template +__global__ void fill_DR_kernel(int m, + const int* __restrict__ csr_row_ptr, + const T* __restrict__ csr_val, + T* __restrict__ DR) +{ + int tid = threadIdx.x; + int bid = blockIdx.x; + int gid = tid + BLOCKSIZE * bid; + + int lid = tid & WARPSIZE - 1; + int wid = tid / WARPSIZE; + + int row = gid / WARPSIZE; + + if(row < m) + { + const int row_start = csr_row_ptr[row]; + const int row_end = csr_row_ptr[row + 1]; + + T max = static_cast(0); + + for(int j = row_start + lid; j < row_end; j += WARPSIZE) + { + const T val = csr_val[j]; + + max = linalg::max(max, linalg::abs(val)); + } + + warp_reduction_max(&max); + + if(lid == 0) + { + DR[row] = max; + } + } +} + + + + +template +__global__ void compute_max_divergence_part1(int m, T eps, T* __restrict__ DR, T* workspace) +{ + int tid = threadIdx.x; + int bid = blockIdx.x; + int gid = tid + BLOCKSIZE * bid; + + __shared__ T shared[BLOCKSIZE]; + + int index = gid; + + T max_divergence = static_cast(0); + + while(index < m) + { + max_divergence = linalg::max(max_divergence, linalg::abs(static_cast(1) - DR[index])); + + DR[index] = (DR[index] > eps) ? std::sqrt(DR[index]) : static_cast(1); + + index += BLOCKSIZE * gridDim.x; + } + + shared[tid] = max_divergence; + __syncthreads(); + + block_reduction_max(&shared[0], tid); + + if(tid == 0) + { + workspace[bid] = shared[0]; + } +} + +template +__global__ void compute_max_divergence_part2(int m, const T* __restrict__ workspace, T* max_divergence) +{ + int tid = threadIdx.x; + int bid = blockIdx.x; + + __shared__ T shared[BLOCKSIZE]; + + shared[tid] = workspace[tid]; + __syncthreads(); + + block_reduction_max(&shared[0], tid); + + if(tid == 0) + { + *max_divergence = shared[0]; + } +} + + + + + +template +__global__ void update_A(int m, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + T* __restrict__ csr_val, + T* __restrict__ DR) +{ + int tid = threadIdx.x; + int bid = blockIdx.x; + int gid = tid + BLOCKSIZE * bid; + + int lid = tid & WARPSIZE - 1; + int wid = tid / WARPSIZE; + + int row = gid / WARPSIZE; + + if(row < m) + { + const int row_start = csr_row_ptr[row]; + const int row_end = csr_row_ptr[row + 1]; + + const T DR_row_val = DR[row]; + + for(int j = row_start + lid; j < row_end; j += WARPSIZE) + { + csr_val[j] /= (DR_row_val * DR[csr_col_ind[j]]); + } + } +} + +template +__global__ void update_D(int m, const T* __restrict__ DR, T* __restrict__ D) +{ + int tid = threadIdx.x; + int bid = blockIdx.x; + int gid = tid + BLOCKSIZE * bid; + + if(gid < m) + { + D[gid] = D[gid] / DR[gid]; + } +} + +#endif diff --git a/library/src/backend/device/device_ruiz_scaling.cpp b/library/src/backend/device/device_ruiz_scaling.cpp index d028179..04d64b3 100644 --- a/library/src/backend/device/device_ruiz_scaling.cpp +++ b/library/src/backend/device/device_ruiz_scaling.cpp @@ -35,13 +35,35 @@ #include "cuda/cuda_ruiz_scaling.h" #endif -void linalg::device_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2) +void linalg::device_ruiz_scaling( + vector& D1, csr_matrix& A, vector& D2, int max_k, double tol) { ROUTINE_TRACE("linalg::device_ruiz_scaling"); if constexpr(is_cuda_available()) { - CALL_CUDA(cuda_ruiz_scaling( - D1.get_vec(), A.get_row_ptr(), A.get_col_ind(), A.get_val(), A.get_m(), D2.get_vec())); + CALL_CUDA(cuda_ruiz_scaling(D1.get_vec(), + A.get_row_ptr(), + A.get_col_ind(), + A.get_val(), + A.get_m(), + D2.get_vec(), + max_k, + tol)); + } + else + { + std::cout << "Error: Not device backend available for the function " << __func__ + << std::endl; + } +} + +void linalg::device_symmetric_ruiz_scaling(vector& D, csr_matrix& A, int max_k, double tol) +{ + ROUTINE_TRACE("linalg::device_symmetric_ruiz_scaling"); + if constexpr(is_cuda_available()) + { + CALL_CUDA(cuda_symmetric_ruiz_scaling( + D.get_vec(), A.get_row_ptr(), A.get_col_ind(), A.get_val(), A.get_m(), max_k, tol)); } else { diff --git a/library/src/backend/device/device_ruiz_scaling.h b/library/src/backend/device/device_ruiz_scaling.h index 1d1381a..64eed10 100644 --- a/library/src/backend/device/device_ruiz_scaling.h +++ b/library/src/backend/device/device_ruiz_scaling.h @@ -31,7 +31,10 @@ namespace linalg { - void device_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2); + void device_ruiz_scaling( + vector& D1, csr_matrix& A, vector& D2, int max_k, double tol); + + void device_symmetric_ruiz_scaling(vector& D, csr_matrix& A, int max_k, double tol); } #endif diff --git a/library/src/backend/host/host_ruiz_scaling.cpp b/library/src/backend/host/host_ruiz_scaling.cpp index 3d21d94..7e3edaa 100644 --- a/library/src/backend/host/host_ruiz_scaling.cpp +++ b/library/src/backend/host/host_ruiz_scaling.cpp @@ -35,25 +35,41 @@ namespace linalg { template - static void host_ruiz_scaling_impl( - T* D1, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, T* D2) + static void host_ruiz_scaling_impl(T* D1, + const int* csr_row_ptr, + const int* csr_col_ind, + T* csr_val, + int m, + T* D2, + int max_k, + T tol) { - ROUTINE_TRACE("host_scale_diagonal_impl"); + ROUTINE_TRACE("host_ruiz_scaling_impl"); - int k = 10; - - // D_1^0 = I, D_2^0 = I + // D1^0 = I, D2^0 = I +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif for(int i = 0; i < m; i++) { D1[i] = static_cast(1); D2[i] = static_cast(1); } - std::vector DR(m, static_cast(0)); - std::vector DC(m, static_cast(0)); + std::vector DR(m); + std::vector DC(m); - for(int iter = 0; iter < k; iter++) + for(int iter = 0; iter < max_k; iter++) { +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif + for(int i = 0; i < m; i++) + { + DR[i] = static_cast(0); + DC[i] = static_cast(0); + } + for(int i = 0; i < m; i++) { const int start = csr_row_ptr[i]; @@ -69,20 +85,29 @@ namespace linalg } } - T row_max = static_cast(0); - T col_max = static_cast(0); + const T eps = std::numeric_limits::epsilon(); + + T max_divergence = static_cast(0); for(int i = 0; i < m; i++) { - row_max = std::max(row_max, DR[i]); - col_max = std::max(col_max, DC[i]); + max_divergence = std::max(max_divergence, std::abs(static_cast(1.0) - DR[i])); + max_divergence = std::max(max_divergence, std::abs(static_cast(1.0) - DC[i])); - DR[i] = std::sqrt(DR[i]); - DC[i] = std::sqrt(DC[i]); + DR[i] = (DR[i] > eps) ? std::sqrt(DR[i]) : static_cast(1); + DC[i] = (DC[i] > eps) ? std::sqrt(DC[i]) : static_cast(1); } - std::cout << "row_max: " << row_max << " col_max: " << col_max << std::endl; + std::cout << "max_divergence: " << max_divergence << std::endl; + + if(max_divergence < tol) + { + break; + } // Update A_k+1 = DR^-1 * A_k * DC^-1 +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif for(int i = 0; i < m; i++) { const int start = csr_row_ptr[i]; @@ -96,6 +121,9 @@ namespace linalg } // Update D1_k+1 = D1_k * DR^-1, D2_k+1 = D2_k * DC^-1 +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif for(int i = 0; i < m; i++) { D1[i] = D1[i] / DR[i]; @@ -103,9 +131,98 @@ namespace linalg } } } + + template + static void host_symmetric_ruiz_scaling_impl( + T* D, const int* csr_row_ptr, const int* csr_col_ind, T* csr_val, int m, int max_k, T tol) + { + ROUTINE_TRACE("host_symmetric_ruiz_scaling_impl"); + + // D^0 = I +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif + for(int i = 0; i < m; i++) + { + D[i] = static_cast(1); + } + + // DR = DC + std::vector DR(m); + + for(int iter = 0; iter < max_k; iter++) + { +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif + for(int i = 0; i < m; i++) + { + DR[i] = static_cast(0); + } + +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif + for(int i = 0; i < m; i++) + { + const int start = csr_row_ptr[i]; + const int end = csr_row_ptr[i + 1]; + + for(int j = start; j < end; j++) + { + const T val = csr_val[j]; + + DR[i] = std::max(DR[i], std::abs(val)); + } + } + + const T eps = std::numeric_limits::epsilon(); + + T max_divergence = static_cast(0); + for(int i = 0; i < m; i++) + { + max_divergence = std::max(max_divergence, std::abs(static_cast(1) - DR[i])); + + DR[i] = (DR[i] > eps) ? std::sqrt(DR[i]) : static_cast(1); + } + + std::cout << "max_divergence: " << max_divergence << std::endl; + + if(max_divergence < tol) + { + break; + } + + // Update A_k+1 = DR^-1 * A_k * DR^-1 +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif + for(int i = 0; i < m; i++) + { + const int start = csr_row_ptr[i]; + const int end = csr_row_ptr[i + 1]; + + for(int j = start; j < end; j++) + { + const int col = csr_col_ind[j]; + csr_val[j] /= (DR[i] * DR[col]); + } + } + + // Update D_k+1 = D_k * DR^-1 +#if defined(_OPENMP) +#pragma omp parallel for schedule(dynamic, 1024) +#endif + for(int i = 0; i < m; i++) + { + D[i] = D[i] / DR[i]; + } + } + } } -void linalg::host_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2) +void linalg::host_ruiz_scaling( + vector& D1, csr_matrix& A, vector& D2, int max_k, double tol) { ROUTINE_TRACE("linalg::host_ruiz_scaling"); @@ -113,6 +230,23 @@ void linalg::host_ruiz_scaling(vector& D1, csr_matrix& A, vector assert(D1.get_size() == A.get_m()); assert(D2.get_size() == A.get_m()); - host_ruiz_scaling_impl( - D1.get_vec(), A.get_row_ptr(), A.get_col_ind(), A.get_val(), A.get_m(), D2.get_vec()); + host_ruiz_scaling_impl(D1.get_vec(), + A.get_row_ptr(), + A.get_col_ind(), + A.get_val(), + A.get_m(), + D2.get_vec(), + max_k, + tol); +} + +void linalg::host_symmetric_ruiz_scaling(vector& D, csr_matrix& A, int max_k, double tol) +{ + ROUTINE_TRACE("linalg::host_ruiz_scaling"); + + assert(A.get_m() == A.get_n()); + assert(D.get_size() == A.get_m()); + + host_symmetric_ruiz_scaling_impl( + D.get_vec(), A.get_row_ptr(), A.get_col_ind(), A.get_val(), A.get_m(), max_k, tol); } diff --git a/library/src/backend/host/host_ruiz_scaling.h b/library/src/backend/host/host_ruiz_scaling.h index 5c338a5..83ed18a 100644 --- a/library/src/backend/host/host_ruiz_scaling.h +++ b/library/src/backend/host/host_ruiz_scaling.h @@ -32,7 +32,10 @@ namespace linalg { - void host_ruiz_scaling(vector& D1, csr_matrix& A, vector& D2); + void host_ruiz_scaling( + vector& D1, csr_matrix& A, vector& D2, int max_k, double tol); + + void host_symmetric_ruiz_scaling(vector& D, csr_matrix& A, int max_k, double tol); } #endif // HOST_APPLY_RUIZ_SCALING_H diff --git a/library/src/csr_matrix.cpp b/library/src/csr_matrix.cpp index 6f1aa97..7f73225 100644 --- a/library/src/csr_matrix.cpp +++ b/library/src/csr_matrix.cpp @@ -418,7 +418,7 @@ void csr_matrix::transpose(csr_matrix& T) const "linalg::csr_matrix::transpose", host_transpose_matrix, device_transpose_matrix, *this, T); } -void csr_matrix::apply_ruiz_scaling(vector& D1, vector& D2) +void csr_matrix::apply_ruiz_scaling(vector& D1, vector& D2, int max_k, double tol) { ROUTINE_TRACE("csr_matrix::apply_ruiz_scaling"); @@ -427,7 +427,22 @@ void csr_matrix::apply_ruiz_scaling(vector& D1, vector& D2) device_ruiz_scaling, D1, *this, - D2); + D2, + max_k, + tol); +} + +void csr_matrix::apply_symmetric_ruiz_scaling(vector& D, int max_k, double tol) +{ + ROUTINE_TRACE("csr_matrix::apply_symmetric_ruiz_scaling"); + + backend_dispatch("linalg::csr_matrix::apply_symmetric_ruiz_scaling", + host_symmetric_ruiz_scaling, + device_symmetric_ruiz_scaling, + D, + *this, + max_k, + tol); } // Structure to hold triplet (COO) format data diff --git a/library/src/iterative_solvers/amg/amg.cpp b/library/src/iterative_solvers/amg/amg.cpp index 633fb12..d82bf87 100644 --- a/library/src/iterative_solvers/amg/amg.cpp +++ b/library/src/iterative_solvers/amg/amg.cpp @@ -108,13 +108,13 @@ namespace linalg static void apply_smoother(const csr_matrix& A, vector& x, const vector& b, - Smoother smoother) + smoother smoother) { ROUTINE_TRACE("apply_smoother"); switch(smoother) { - case Smoother::Jacobi: + case smoother::jacobi: { vector xold(A.get_m()); xold.copy_from(x); @@ -122,16 +122,16 @@ namespace linalg jacobi_iteration(A, x, xold, b); break; } - case Smoother::Gauss_Seidel: + case smoother::gauss_seidel: gauss_seidel_iteration(A, x, b); break; - case Smoother::Symm_Gauss_Seidel: + case smoother::symmetric_gauss_seidel: gauss_seidel_iteration(A, x, b); break; - case Smoother::SOR: + case smoother::SOR: sor_iteration(A, x, b, 0.5f); break; - case Smoother::SSOR: + case smoother::SSOR: ssor_iteration(A, x, b, 0.5f); break; default: @@ -145,7 +145,7 @@ namespace linalg int n1, int n2, int currentLevel, - Smoother smoother) + smoother smoother) { ROUTINE_TRACE("vcycle"); @@ -271,7 +271,7 @@ namespace linalg int n2, int n3, int currentLevel, - Smoother smoother) + smoother smoother) { ROUTINE_TRACE("wcycle"); @@ -394,7 +394,7 @@ namespace linalg int n2, int n3, int currentLevel, - Smoother smoother) + smoother smoother) { ROUTINE_TRACE("fcycle"); @@ -516,8 +516,8 @@ int linalg::amg_solve(const hierarchy& hierarchy, const vector& b, int n1, int n2, - Cycle cycle, - Smoother smoother, + cycle cycle, + smoother smoother, iter_control control) { ROUTINE_TRACE("amg_solve"); @@ -538,13 +538,13 @@ int linalg::amg_solve(const hierarchy& hierarchy, { switch(cycle) { - case Cycle::Vcycle: + case cycle::vcycle: vcycle(hierarchy, x, b, n1, n2, 0, smoother); break; - case Cycle::Wcycle: + case cycle::wcycle: wcycle(hierarchy, x, b, n1, n2, n3, 0, smoother); break; - case Cycle::Fcycle: + case cycle::fcycle: fcycle(hierarchy, x, b, n1, n2, n3, 0, smoother); break; } @@ -566,4 +566,4 @@ int linalg::amg_solve(const hierarchy& hierarchy, std::cout << "amg solve time: " << ms_double.count() << "ms" << std::endl; return cycle_count; -} \ No newline at end of file +} From ca147ea38b766554cd423476e3aa2ac1dba90ebf Mon Sep 17 00:00:00 2001 From: jsandham Date: Thu, 30 Jul 2026 18:47:49 -0400 Subject: [PATCH 3/3] Add missing limits and cmath includes --- library/src/backend/host/host_ruiz_scaling.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/library/src/backend/host/host_ruiz_scaling.cpp b/library/src/backend/host/host_ruiz_scaling.cpp index 7e3edaa..8265bcc 100644 --- a/library/src/backend/host/host_ruiz_scaling.cpp +++ b/library/src/backend/host/host_ruiz_scaling.cpp @@ -26,7 +26,9 @@ #include #include +#include #include +#include #include "host_ruiz_scaling.h"