diff --git a/common/cuda_hip/CMakeLists.txt b/common/cuda_hip/CMakeLists.txt index 3ecf62951f0..69a9a84c25a 100644 --- a/common/cuda_hip/CMakeLists.txt +++ b/common/cuda_hip/CMakeLists.txt @@ -36,6 +36,7 @@ set(CUDA_HIP_SOURCES matrix/sellp_kernels.cpp matrix/sparsity_csr_kernels.cpp multigrid/pgm_kernels.cpp + multigrid/pmis_kernels.cpp preconditioner/isai_kernels.cpp preconditioner/jacobi_kernels.cpp preconditioner/jacobi_advanced_apply_kernels.cpp diff --git a/common/cuda_hip/multigrid/pmis_kernels.cpp b/common/cuda_hip/multigrid/pmis_kernels.cpp new file mode 100644 index 00000000000..d6d788308f7 --- /dev/null +++ b/common/cuda_hip/multigrid/pmis_kernels.cpp @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/multigrid/pmis_kernels.hpp" + +#include + +#include + +#include "common/cuda_hip/base/randlib_bindings.hpp" + +namespace gko { +namespace kernels { +namespace GKO_DEVICE_NAMESPACE { +namespace pmis { + + +template +void initialize_random_weight(std::shared_ptr exec, + size_type num, ValueType* weight) +{ + auto gen = randlib::rand_generator( + std::random_device{}(), RANDLIB_RNG_PSEUDO_DEFAULT, exec->get_stream()); + randlib::uniform_rand_vector(gen, num, weight); + randlib::destroy(gen); +} + +GKO_INSTANTIATE_FOR_EACH_NON_COMPLEX_VALUE_TYPE_BASE( + GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL); + + +} // namespace pmis +} // namespace GKO_DEVICE_NAMESPACE +} // namespace kernels +} // namespace gko diff --git a/common/unified/multigrid/pmis_kernels.cpp b/common/unified/multigrid/pmis_kernels.cpp index 27af79540b8..cbce3be2495 100644 --- a/common/unified/multigrid/pmis_kernels.cpp +++ b/common/unified/multigrid/pmis_kernels.cpp @@ -27,12 +27,32 @@ namespace GKO_DEVICE_NAMESPACE { namespace pmis { +// the number of threads working on the same row +constexpr int width = 32; + + template void compute_row_maxabs(std::shared_ptr exec, const matrix::Csr* csr, remove_complex* row_maxabs) { - GKO_NOT_IMPLEMENTED; + run_kernel_row_reduction( + exec, + [] GKO_KERNEL(auto row, auto tid, auto row_ptrs, auto col_idxs, + auto values) { + auto maxabs = zero(abs(values[0])); + for (auto idx = tid + row_ptrs[row]; idx < row_ptrs[row + 1]; + idx += width) { + if (row == col_idxs[idx]) { + continue; + } + maxabs = gko::max(maxabs, abs(values[idx])); + } + return maxabs; + }, + GKO_KERNEL_REDUCE_MAX(remove_complex), row_maxabs, 1, + dim<2>{csr->get_size()[0], width}, csr->get_const_row_ptrs(), + csr->get_const_col_idxs(), csr->get_const_values()); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( @@ -46,7 +66,31 @@ void compute_strong_dep_row(std::shared_ptr exec, remove_complex strength_threshold, IndexType* sparsity_rows) { - GKO_NOT_IMPLEMENTED; + run_kernel_row_reduction( + exec, + [] GKO_KERNEL(auto row, auto tid, auto row_maxabs, + auto strength_threshold, auto row_ptrs, auto col_idxs, + auto values) { + auto max_abs = row_maxabs[row]; + auto count = zero(); + if (max_abs == zero(max_abs)) { + return count; + } + for (auto idx = tid + row_ptrs[row]; idx < row_ptrs[row + 1]; + idx += width) { + if (row == col_idxs[idx]) { + continue; + } + if (abs(values[idx]) >= strength_threshold * max_abs) { + count++; + } + } + return count; + }, + GKO_KERNEL_REDUCE_SUM(IndexType), sparsity_rows, 1, + dim<2>{csr->get_size()[0], width}, row_maxabs, strength_threshold, + csr->get_const_row_ptrs(), csr->get_const_col_idxs(), + csr->get_const_values()); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( @@ -60,7 +104,33 @@ void compute_strong_dep(std::shared_ptr exec, remove_complex strength_threshold, matrix::SparsityCsr* strong_dep) { - GKO_NOT_IMPLEMENTED; + // we handle this by one thread per row. It might get improved if we use a + // warp with popcount and prefix for a row. + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto row_maxabs, auto strength_threshold, + auto row_ptrs, auto col_idxs, auto values, + auto dep_row_ptrs, auto dep_col_idxs) { + auto max_abs = row_maxabs[row]; + if (max_abs == zero(max_abs)) { + return; + } + auto d_idx = dep_row_ptrs[row]; + for (auto idx = row_ptrs[row]; idx < row_ptrs[row + 1]; idx++) { + const auto col = col_idxs[idx]; + if (row == col) { + continue; + } + if (abs(values[idx]) >= strength_threshold * max_abs) { + dep_col_idxs[d_idx] = col; + d_idx++; + } + } + }, + csr->get_size()[0], row_maxabs, strength_threshold, + csr->get_const_row_ptrs(), csr->get_const_col_idxs(), + csr->get_const_values(), strong_dep->get_const_row_ptrs(), + strong_dep->get_col_idxs()); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( @@ -73,7 +143,30 @@ void initialize_weight_and_status( const matrix::SparsityCsr* trans_strong_dep, remove_complex* weight, int* status) { - GKO_NOT_IMPLEMENTED; + auto num = trans_strong_dep->get_size()[0]; + array random(exec, num); + // note. range setting 0, 1 has different meaning in different backend + // for range(l, r) + // std include `l` but exclude `r` + // cuda/hip exclude `l` but include `r` + // dpcpp does not mentioned it in the documentation but the code should + // include `l` but exclude `r`. We only use float here because dpcpp device + // may lack of double precision support and random generator does not + // support 16-bit. + initialize_random_weight(exec, num, random.get_data()); + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto row_ptrs, auto random, auto weight, + auto status) { + using type = device_type>; + auto w = static_cast(row_ptrs[row + 1] - row_ptrs[row]); + status[row] = + (w == 0.0f ? kernels::pmis::fine : kernels::pmis::unassigned); + // avoid random value to be 1 + weight[row] = static_cast(random[row] * 0.99f + w); + }, + num, trans_strong_dep->get_const_row_ptrs(), random.get_const_data(), + weight, status); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( @@ -84,10 +177,60 @@ template void classify(std::shared_ptr exec, const remove_complex* weight, const matrix::SparsityCsr* strong_dep, - const matrix::SparsityCsr* trans_strong_dep, const int* status, int* new_status) { - GKO_NOT_IMPLEMENTED; + static_assert(kernels::pmis::unassigned < kernels::pmis::coarse, + "we use min reduction to mark local maximum as coarse"); + // mark coarse point + run_kernel_row_reduction( + exec, + [] GKO_KERNEL(auto row, auto tid, auto status, auto weight, + auto row_ptrs, auto col_idxs) { + auto ans = status[row]; + if (ans != kernels::pmis::unassigned) { + return ans; + } + for (auto idx = tid + row_ptrs[row]; idx < row_ptrs[row + 1]; + idx += width) { + auto col = col_idxs[idx]; + if (status[col] == kernels::pmis::unassigned && + device_std::tie(weight[col], col) > + device_std::tie(weight[row], row)) { + return kernels::pmis::unassigned; + } + } + return kernels::pmis::coarse; + }, + [] GKO_KERNEL(auto a, auto b) { return a < b ? a : b; } /* minimun */, + [] GKO_KERNEL(auto a) { return a; }, int{1}, new_status, 1, + dim<2>{strong_dep->get_size()[0], width}, status, weight, + strong_dep->get_const_row_ptrs(), strong_dep->get_const_col_idxs()); + // mark new fine point strongly influenced by the new coarse points + // TODO: using warp vote function if implement in native way. + static_assert(kernels::pmis::fine > kernels::pmis::unassigned, + "we use max reduction to mark new fine by any strong coarse"); + run_kernel_row_reduction( + exec, + [] GKO_KERNEL(auto row, auto tid, auto new_status, auto row_ptrs, + auto col_idxs) { + if (new_status[row] != kernels::pmis::unassigned) { + return new_status[row]; + } + for (auto idx = tid + row_ptrs[row]; idx < row_ptrs[row + 1]; + idx += width) { + // we will only update new_status from -1 to 0 or keep -1, so + // grabbing this value is fine no matter if it is updated or + // not. + if (new_status[col_idxs[idx]] == kernels::pmis::coarse) { + return kernels::pmis::fine; + } + } + return kernels::pmis::unassigned; + }, + [] GKO_KERNEL(auto a, auto b) { return a > b ? a : b; } /* maximum */, + [] GKO_KERNEL(auto a) { return a; }, int{-1}, new_status, 1, + dim<2>{strong_dep->get_size()[0], width}, new_status, + strong_dep->get_const_row_ptrs(), strong_dep->get_const_col_idxs()); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_PMIS_CLASSIFY_KERNEL); @@ -96,7 +239,15 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_PMIS_CLASSIFY_KERNEL); void count(std::shared_ptr exec, size_type num, const int* status, size_type* num_unassigned) { - GKO_NOT_IMPLEMENTED; + array d_result(exec, 1); + run_kernel_reduction( + exec, + [] GKO_KERNEL(auto i, auto status) { + return static_cast(status[i] == + kernels::pmis::unassigned); + }, + GKO_KERNEL_REDUCE_SUM(size_type), d_result.get_data(), num, status); + *num_unassigned = get_element(d_result, 0); } @@ -106,7 +257,25 @@ void direct_interpolation_row_count( const matrix::SparsityCsr* strong_dep, const int* status, IndexType* prolong_row_ptr) { - GKO_NOT_IMPLEMENTED; + run_kernel_row_reduction( + exec, + [] GKO_KERNEL(auto row, auto tid, auto status, auto row_ptrs, + auto col_idxs) { + if (status[row] == kernels::pmis::coarse) { + return tid == 0 ? one() : zero(); + } + auto count = zero(); + for (auto idx = tid + row_ptrs[row]; idx < row_ptrs[row + 1]; + idx += width) { + if (status[col_idxs[idx]] == kernels::pmis::coarse) { + count++; + } + } + return count; + }, + GKO_KERNEL_REDUCE_SUM(IndexType), prolong_row_ptr, 1, + dim<2>{strong_dep->get_size()[0], width}, status, + strong_dep->get_const_row_ptrs(), strong_dep->get_const_col_idxs()); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( @@ -122,7 +291,83 @@ void direct_interpolation_fill( const IndexType* coarse_map, const IndexType* prolong_row_ptrs, IndexType* prolong_col_idxs, ValueType* prolong_values) { - GKO_NOT_IMPLEMENTED; + // currently use one thread per row. It might get improved by using a warp + // for row with prefix and popcount + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto row_maxabs, auto strength_threshold, + auto coarse_map, auto row_ptrs, auto col_idxs, + auto values, auto prolong_row_ptrs, auto prolong_col_idxs, + auto prolong_values) { + if (coarse_map[row] != coarse_map[row + 1]) { + auto idx = prolong_row_ptrs[row]; + prolong_col_idxs[idx] = coarse_map[row]; + prolong_values[idx] = one(prolong_values[idx]); + return; + } + auto pos = zero(values[0]); + auto pos_divisor = zero(values[0]); + auto neg = zero(values[0]); + auto neg_divisor = zero(values[0]); + auto diag = zero(values[0]); + bool enable_neg = false; + bool enable_pos = false; + // first compute alpha/beta + auto max_abs = row_maxabs[row]; + for (auto idx = row_ptrs[row]; idx < row_ptrs[row + 1]; idx++) { + auto val = values[idx]; + auto col = col_idxs[idx]; + if (col == row) { + diag = val; + continue; + } + if (real(val) >= zero(real(val))) { + pos += val; + if (coarse_map[col] != coarse_map[col + 1] && + abs(val) >= strength_threshold * max_abs) { + pos_divisor += val; + enable_pos = true; + } + } else { + neg += val; + if (coarse_map[col] != coarse_map[col + 1] && + abs(val) >= strength_threshold * max_abs) { + neg_divisor += val; + enable_neg = true; + } + } + } + pos = safe_divide(pos, pos_divisor); + neg = safe_divide(neg, neg_divisor); + if (!enable_neg && !enable_pos) { + return; + } + + auto p_idx = prolong_row_ptrs[row]; + for (auto idx = row_ptrs[row]; idx < row_ptrs[row + 1]; idx++) { + auto val = values[idx]; + auto col = col_idxs[idx]; + if (col == row || abs(val) < strength_threshold * max_abs) { + continue; + } + if (real(val) >= zero(real(val)) && enable_pos && + coarse_map[col] != coarse_map[col + 1]) { + prolong_col_idxs[p_idx] = coarse_map[col]; + prolong_values[p_idx] = -pos * val / diag; + p_idx++; + } + if (real(val) < zero(real(val)) && enable_neg && + coarse_map[col] != coarse_map[col + 1]) { + prolong_col_idxs[p_idx] = coarse_map[col]; + prolong_values[p_idx] = -neg * val / diag; + p_idx++; + } + } + }, + csr->get_size()[0], row_maxabs, strength_threshold, coarse_map, + csr->get_const_row_ptrs(), csr->get_const_col_idxs(), + csr->get_const_values(), prolong_row_ptrs, prolong_col_idxs, + prolong_values); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( diff --git a/core/device_hooks/common_kernels.inc.cpp b/core/device_hooks/common_kernels.inc.cpp index 28bd55534d4..5286f197800 100644 --- a/core/device_hooks/common_kernels.inc.cpp +++ b/core/device_hooks/common_kernels.inc.cpp @@ -1153,6 +1153,8 @@ namespace pmis { GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_PMIS_COMPUTE_ROW_MAXABS_KERNEL); GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_PMIS_COMPUTE_STRONG_DEP_ROW_KERNEL); GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_PMIS_COMPUTE_STRONG_DEP_KERNEL); +GKO_STUB_NON_COMPLEX_VALUE_TYPE_BASE( + GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL); GKO_STUB_VALUE_AND_INDEX_TYPE( GKO_DECLARE_PMIS_INITIALIZE_WEIGHT_AND_STATUS_KERNEL); GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_PMIS_CLASSIFY_KERNEL); diff --git a/core/multigrid/pmis.cpp b/core/multigrid/pmis.cpp index a662caf4e18..a24b6b95f99 100644 --- a/core/multigrid/pmis.cpp +++ b/core/multigrid/pmis.cpp @@ -163,9 +163,9 @@ void Pmis::generate() exec->run( pmis::make_count(this->get_size()[0], status_ptr, &num_not_assigned)); while (num_not_assigned != 0) { - exec->run(pmis::make_classify( - weight_.get_const_data(), strong_dep.get(), - transpose_strong_dep.get(), status_ptr, new_status_ptr)); + exec->run(pmis::make_classify(weight_.get_const_data(), + strong_dep.get(), status_ptr, + new_status_ptr)); size_type new_num = 0; exec->run( pmis::make_count(this->get_size()[0], new_status_ptr, &new_num)); diff --git a/core/multigrid/pmis_kernels.hpp b/core/multigrid/pmis_kernels.hpp index cdb68c4b37f..d3253a86de4 100644 --- a/core/multigrid/pmis_kernels.hpp +++ b/core/multigrid/pmis_kernels.hpp @@ -46,6 +46,10 @@ constexpr int unassigned = -1; remove_complex strength_threshold, \ matrix::SparsityCsr* strong_dep) +#define GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL(ValueType) \ + void initialize_random_weight(std::shared_ptr exec, \ + size_type num, ValueType* weight) + #define GKO_DECLARE_PMIS_INITIALIZE_WEIGHT_AND_STATUS_KERNEL(ValueType, \ IndexType) \ void initialize_weight_and_status( \ @@ -53,13 +57,11 @@ constexpr int unassigned = -1; const matrix::SparsityCsr* trans_strong_dep, \ remove_complex* weight, int* status) -#define GKO_DECLARE_PMIS_CLASSIFY_KERNEL(ValueType, IndexType) \ - void classify( \ - std::shared_ptr exec, \ - const remove_complex* weight, \ - const matrix::SparsityCsr* strong_dep, \ - const matrix::SparsityCsr* trans_strong_dep, \ - const int* status, int* new_status) +#define GKO_DECLARE_PMIS_CLASSIFY_KERNEL(ValueType, IndexType) \ + void classify(std::shared_ptr exec, \ + const remove_complex* weight, \ + const matrix::SparsityCsr* strong_dep, \ + const int* status, int* new_status) #define GKO_DECLARE_COUNT_KERNEL \ void count(std::shared_ptr exec, size_type num, \ @@ -88,6 +90,8 @@ constexpr int unassigned = -1; GKO_DECLARE_PMIS_COMPUTE_STRONG_DEP_ROW_KERNEL(ValueType, IndexType); \ template \ GKO_DECLARE_PMIS_COMPUTE_STRONG_DEP_KERNEL(ValueType, IndexType); \ + template \ + GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL(ValueType); \ template \ GKO_DECLARE_PMIS_INITIALIZE_WEIGHT_AND_STATUS_KERNEL(ValueType, \ IndexType); \ diff --git a/cuda/base/curand_bindings.hpp b/cuda/base/curand_bindings.hpp index 80ceff2dacd..f1e1eeab241 100644 --- a/cuda/base/curand_bindings.hpp +++ b/cuda/base/curand_bindings.hpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2026 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -95,6 +95,25 @@ GKO_BIND_CURAND_RANDOM_VECTOR(ValueType, detail::not_implemented); #undef GKO_BIND_CURAND_RANDOM_VECTOR +#define GKO_BIND_CURAND_UNIFORM_RANDOM_VECTOR(ValueType, CurandName) \ + inline void uniform_rand_vector(curandGenerator_t& gen, size_type n, \ + ValueType* values) \ + { \ + GKO_ASSERT_NO_CURAND_ERRORS(CurandName(gen, values, n)); \ + } \ + static_assert(true, \ + "This assert is used to counter the false positive extra " \ + "semi-colon warnings") + +GKO_BIND_CURAND_UNIFORM_RANDOM_VECTOR(float, curandGenerateUniform); +GKO_BIND_CURAND_UNIFORM_RANDOM_VECTOR(double, curandGenerateUniformDouble); +template +GKO_BIND_CURAND_UNIFORM_RANDOM_VECTOR(ValueType, detail::not_implemented); + + +#undef GKO_BIND_CURAND_UNIFORM_RANDOM_VECTOR + + } // namespace curand diff --git a/dpcpp/CMakeLists.txt b/dpcpp/CMakeLists.txt index 1abb09c39a6..777271ef594 100644 --- a/dpcpp/CMakeLists.txt +++ b/dpcpp/CMakeLists.txt @@ -72,6 +72,7 @@ target_sources( matrix/sellp_kernels.dp.cpp matrix/sparsity_csr_kernels.dp.cpp multigrid/pgm_kernels.dp.cpp + multigrid/pmis_kernels.dp.cpp preconditioner/batch_jacobi_kernels.dp.cpp preconditioner/isai_kernels.dp.cpp preconditioner/jacobi_advanced_apply_kernel.dp.cpp diff --git a/dpcpp/multigrid/pmis_kernels.dp.cpp b/dpcpp/multigrid/pmis_kernels.dp.cpp new file mode 100644 index 00000000000..c4725062144 --- /dev/null +++ b/dpcpp/multigrid/pmis_kernels.dp.cpp @@ -0,0 +1,46 @@ +// SPDX-FileCopyrightText: 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include "core/multigrid/pmis_kernels.hpp" + +#include + +#include + +#include + +#include "common/cuda_hip/base/randlib_bindings.hpp" + +namespace gko { +namespace kernels { +namespace GKO_DEVICE_NAMESPACE { +namespace pmis { + + +template +void initialize_random_weight(std::shared_ptr exec, + size_type num, ValueType* weight) +{ + auto seed = std::random_device{}(); + exec->get_queue()->submit([&](sycl::handler& cgh) { + cgh.parallel_for(sycl::range<1>(num), [=](sycl::item<1> idx) { + std::uint64_t offset = idx.get_linear_id(); + oneapi::dpl::minstd_rand engine(seed, offset); + oneapi::dpl::uniform_real_distribution> + distr(0, 1); + work[idx] = distr(engine); + }); + }); +} + +GKO_INSTANTIATE_FOR_EACH_NON_COMPLEX_VALUE_TYPE_BASE( + GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL); + + +} // namespace pmis +} // namespace GKO_DEVICE_NAMESPACE +} // namespace kernels +} // namespace gko diff --git a/hip/base/hiprand_bindings.hip.hpp b/hip/base/hiprand_bindings.hip.hpp index f4f4d45f65b..44624330443 100644 --- a/hip/base/hiprand_bindings.hip.hpp +++ b/hip/base/hiprand_bindings.hip.hpp @@ -96,6 +96,25 @@ GKO_BIND_HIPRAND_RANDOM_VECTOR(ValueType, detail::not_implemented); #undef GKO_BIND_HIPRAND_RANDOM_VECTOR +#define GKO_BIND_HIPRAND_UNIFORM_RANDOM_VECTOR(ValueType, HiprandName) \ + inline void uniform_rand_vector(hiprandGenerator_t& gen, size_type n, \ + ValueType* values) \ + { \ + GKO_ASSERT_NO_HIPRAND_ERRORS(HiprandName(gen, values, n)); \ + } \ + static_assert(true, \ + "This assert is used to counter the false positive extra " \ + "semi-colon warnings") + +GKO_BIND_HIPRAND_UNIFORM_RANDOM_VECTOR(float, hiprandGenerateUniform); +GKO_BIND_HIPRAND_UNIFORM_RANDOM_VECTOR(double, hiprandGenerateUniformDouble); +template +GKO_BIND_HIPRAND_UNIFORM_RANDOM_VECTOR(ValueType, detail::not_implemented); + + +#undef GKO_BIND_HIPRAND_UNIFORM_RANDOM_VECTOR + + } // namespace hiprand diff --git a/omp/CMakeLists.txt b/omp/CMakeLists.txt index b7b5c209342..161e899ed74 100644 --- a/omp/CMakeLists.txt +++ b/omp/CMakeLists.txt @@ -46,6 +46,7 @@ target_sources( matrix/sellp_kernels.cpp matrix/sparsity_csr_kernels.cpp multigrid/pgm_kernels.cpp + multigrid/pmis_kernels.cpp preconditioner/batch_jacobi_kernels.cpp preconditioner/isai_kernels.cpp preconditioner/jacobi_kernels.cpp diff --git a/omp/multigrid/pmis_kernels.cpp b/omp/multigrid/pmis_kernels.cpp new file mode 100644 index 00000000000..22a6e6192af --- /dev/null +++ b/omp/multigrid/pmis_kernels.cpp @@ -0,0 +1,34 @@ +// SPDX-FileCopyrightText: 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/multigrid/pmis_kernels.hpp" + +#include + +#include + +namespace gko { +namespace kernels { +namespace omp { +namespace pmis { + + +template +void initialize_random_weight(std::shared_ptr exec, + size_type num, ValueType* weight) +{ + std::default_random_engine gen(42); + std::uniform_real_distribution dist(0.0, 1.0); + for (size_type row = 0; row < num; row++) { + weight[row] = dist(gen); + } +} +GKO_INSTANTIATE_FOR_EACH_NON_COMPLEX_VALUE_TYPE_BASE( + GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL); + + +} // namespace pmis +} // namespace omp +} // namespace kernels +} // namespace gko diff --git a/reference/multigrid/pmis_kernels.cpp b/reference/multigrid/pmis_kernels.cpp index aa893f973d0..0e46cab1930 100644 --- a/reference/multigrid/pmis_kernels.cpp +++ b/reference/multigrid/pmis_kernels.cpp @@ -122,6 +122,21 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( GKO_DECLARE_PMIS_COMPUTE_STRONG_DEP_KERNEL); +template +void initialize_random_weight(std::shared_ptr exec, + size_type num, ValueType* weight) +{ + std::default_random_engine gen(42); + std::uniform_real_distribution dist(0.0, 1.0); + for (size_type row = 0; row < num; row++) { + weight[row] = dist(gen); + } +} + +GKO_INSTANTIATE_FOR_EACH_NON_COMPLEX_VALUE_TYPE_BASE( + GKO_DECLARE_PMIS_INITIALIZE_RANDOM_WEIGHT_KERNEL); + + template void initialize_weight_and_status( std::shared_ptr exec, @@ -153,7 +168,6 @@ template void classify(std::shared_ptr exec, const remove_complex* weight, const matrix::SparsityCsr* strong_dep, - const matrix::SparsityCsr* trans_strong_dep, const int* status, int* new_status) { const auto nrows = static_cast(strong_dep->get_size()[0]); @@ -182,18 +196,15 @@ void classify(std::shared_ptr exec, } new_status[row] = ans; } - // mark all points strongly influenced by the new coarse points to fine - // group + // mark new fine point strongly influenced by the new coarse points for (IndexType row = 0; row < nrows; row++) { - if (new_status[row] == kernels::pmis::coarse && - new_status[row] != status[row]) { - for (auto idx = trans_strong_dep->get_const_row_ptrs()[row]; - idx < trans_strong_dep->get_const_row_ptrs()[row + 1]; idx++) { - // It is correct even if more than one threads might assign the - // value - auto col = trans_strong_dep->get_const_col_idxs()[idx]; - if (new_status[col] == kernels::pmis::unassigned) { - new_status[col] = kernels::pmis::fine; + if (new_status[row] == kernels::pmis::unassigned) { + for (auto idx = strong_dep->get_const_row_ptrs()[row]; + idx < strong_dep->get_const_row_ptrs()[row + 1]; idx++) { + if (new_status[strong_dep->get_const_col_idxs()[idx]] == + kernels::pmis::coarse) { + new_status[row] = kernels::pmis::fine; + break; } } } diff --git a/reference/test/multigrid/pmis_kernels.cpp b/reference/test/multigrid/pmis_kernels.cpp index 6f945edb1cf..c9c0ef1f40d 100644 --- a/reference/test/multigrid/pmis_kernels.cpp +++ b/reference/test/multigrid/pmis_kernels.cpp @@ -239,15 +239,13 @@ TYPED_TEST(Pmis, Classify) SparsityCsr::create(this->exec, this->mtx.at(i)->get_size(), std::move(this->dep_col_idxs.at(i)), std::move(this->dep_row_ptrs.at(i))); - auto trans_strong_dep = gko::as(strong_dep->transpose()); auto new_status = this->expected_status.at(i); for (int step = 0; step < required_step.at(i); step++) { SCOPED_TRACE(step); auto status = new_status; gko::kernels::reference::pmis::classify( this->exec, weight.at(i).get_data(), strong_dep.get(), - trans_strong_dep.get(), status.get_const_data(), - new_status.get_data()); + status.get_const_data(), new_status.get_data()); GKO_ASSERT_ARRAY_EQ(new_status, status_ans.at(status_idx)); status_idx++; @@ -266,13 +264,12 @@ TYPED_TEST(Pmis, ClassifyOnSameWeight) SparsityCsr::create(this->exec, this->mtx.at(0)->get_size(), std::move(this->dep_col_idxs.at(0)), std::move(this->dep_row_ptrs.at(0))); - auto trans_strong_dep = gko::as(strong_dep->transpose()); gko::array status_ans(this->exec, {f, c, c, f}); auto new_status = this->expected_status.at(0); auto status = new_status; gko::kernels::reference::pmis::classify( - this->exec, weight.get_data(), strong_dep.get(), trans_strong_dep.get(), + this->exec, weight.get_data(), strong_dep.get(), status.get_const_data(), new_status.get_data()); GKO_ASSERT_ARRAY_EQ(new_status, status_ans); diff --git a/test/multigrid/CMakeLists.txt b/test/multigrid/CMakeLists.txt index 31b373adc52..f19dcc32de4 100644 --- a/test/multigrid/CMakeLists.txt +++ b/test/multigrid/CMakeLists.txt @@ -1,4 +1,5 @@ ginkgo_create_common_test(pgm_kernels) +ginkgo_create_common_test(pmis_kernels) ginkgo_create_common_test(fixed_coarsening_kernels) ginkgo_create_common_test(uniform_coarsening_kernels) ginkgo_create_common_test(rs_kernels) diff --git a/test/multigrid/pmis_kernels.cpp b/test/multigrid/pmis_kernels.cpp new file mode 100644 index 00000000000..6f308ad38a3 --- /dev/null +++ b/test/multigrid/pmis_kernels.cpp @@ -0,0 +1,314 @@ +// SPDX-FileCopyrightText: 2017 - 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/multigrid/pmis_kernels.hpp" + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "core/components/precision_conversion_kernels.hpp" +#include "core/components/prefix_sum_kernels.hpp" +#include "core/test/utils.hpp" +#include "core/test/utils/matrix_generator.hpp" +#include "core/test/utils/unsort_matrix.hpp" +#include "core/utils/matrix_utils.hpp" +#include "test/utils/common_fixture.hpp" + + +class Pmis : public CommonTestFixture { +protected: + using Csr = gko::matrix::Csr; + using SparsityCsr = gko::matrix::SparsityCsr; + using real_type = gko::remove_complex; + + Pmis() + : rand_engine(30), + row_ptrs(ref), + maxabs(ref), + col_idxs(ref), + weight(ref), + status(ref), + new_status(ref), + final_status(ref), + prolong_row_ptrs(ref), + coarse_map(ref) + {} + + void initialize_data() + { +#ifdef GINKGO_FAST_TESTS + m = 129; +#else + m = 597; +#endif + auto system_data = + gko::test::generate_random_matrix_data( + m, m, std::uniform_int_distribution<>(10, m), + std::normal_distribution(-1.0, 1.0), rand_engine); + gko::utils::make_diag_dominant(system_data); + system_mtx = Csr::create(ref); + system_mtx->read(system_data); + d_system_mtx = gko::clone(exec, system_mtx); + // the followings run almost whole pmis on reference. + auto num = system_mtx->get_size()[0]; + row_ptrs.resize_and_reset(num + 1); + maxabs.resize_and_reset(num); + gko::kernels::reference::pmis::compute_row_maxabs(ref, system_mtx.get(), + maxabs.get_data()); + gko::kernels::reference::pmis::compute_strong_dep_row( + ref, system_mtx.get(), maxabs.get_const_data(), real_type{0.25}, + row_ptrs.get_data()); + gko::kernels::reference::components::prefix_sum_nonnegative( + ref, row_ptrs.get_data(), row_ptrs.get_size()); + col_idxs.resize_and_reset(row_ptrs.get_const_data()[num]); + strong_dep = gko::matrix::SparsityCsr::create( + ref, system_mtx->get_size(), col_idxs, row_ptrs); + gko::kernels::reference::pmis::compute_strong_dep( + ref, system_mtx.get(), maxabs.get_const_data(), real_type{0.25}, + strong_dep.get()); + trans_strong_dep = gko::as(strong_dep->transpose()); + weight.resize_and_reset(num); + status.resize_and_reset(num); + gko::kernels::reference::pmis::initialize_weight_and_status( + ref, trans_strong_dep.get(), weight.get_data(), status.get_data()); + new_status.resize_and_reset(num); + auto status_ptr = status.get_data(); + auto new_status_ptr = new_status.get_data(); + gko::size_type num_not_assigned = 0; + gko::kernels::reference::pmis::count(ref, num, status_ptr, + &num_not_assigned); + while (num_not_assigned != 0) { + gko::kernels::reference::pmis::classify(ref, weight.get_data(), + strong_dep.get(), + status_ptr, new_status_ptr); + gko::size_type new_num = 0; + gko::kernels::reference::pmis::count(ref, num, new_status_ptr, + &new_num); + if (new_num == num_not_assigned) { + // no progess -> throw error (maybe unneccessary) + throw std::runtime_error("no progress in Pmis"); + } + num_not_assigned = new_num; + std::swap(new_status_ptr, status_ptr); + } + if (status_ptr == status.get_data()) { + final_status = status; + } else { + final_status = new_status; + } + + prolong_row_ptrs.resize_and_reset(num + 1); + gko::kernels::reference::pmis::direct_interpolation_row_count( + ref, strong_dep.get(), status_ptr, prolong_row_ptrs.get_data()); + gko::kernels::reference::components::prefix_sum_nonnegative( + ref, prolong_row_ptrs.get_data(), prolong_row_ptrs.get_size()); + coarse_map.resize_and_reset(num + 1); + gko::kernels::reference::components::convert_precision( + ref, num, status_ptr, coarse_map.get_data()); + gko::kernels::reference::components::prefix_sum_nonnegative( + ref, coarse_map.get_data(), coarse_map.get_size()); + auto prolong_nnz = prolong_row_ptrs.get_const_data()[num]; + } + + std::default_random_engine rand_engine; + std::shared_ptr system_mtx; + std::shared_ptr d_system_mtx; + gko::size_type m; + gko::array row_ptrs; + gko::array maxabs; + gko::array col_idxs; + std::shared_ptr strong_dep; + std::shared_ptr trans_strong_dep; + gko::array weight; + gko::array status; + gko::array new_status; + gko::array final_status; + gko::array prolong_row_ptrs; + gko::array coarse_map; +}; + + +TEST_F(Pmis, ComputeRowMaxAbsIsEquivalentToRef) +{ + initialize_data(); + gko::array maxabs(ref, system_mtx->get_size()[0]); + gko::array d_maxabs(exec, d_system_mtx->get_size()[0]); + + gko::kernels::reference::pmis::compute_row_maxabs(ref, system_mtx.get(), + maxabs.get_data()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::compute_row_maxabs( + exec, d_system_mtx.get(), d_maxabs.get_data()); + + GKO_ASSERT_ARRAY_NEAR(d_maxabs, maxabs, r::value); +} + + +TEST_F(Pmis, ComputeStrongDepRowIsEquivalentToRef) +{ + initialize_data(); + gko::array d_maxabs(exec, maxabs); + gko::array rows(ref, system_mtx->get_size()[0]); + gko::array d_rows(exec, d_system_mtx->get_size()[0]); + + gko::kernels::reference::pmis::compute_strong_dep_row( + ref, system_mtx.get(), maxabs.get_const_data(), real_type{0.25}, + rows.get_data()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::compute_strong_dep_row( + exec, d_system_mtx.get(), d_maxabs.get_const_data(), real_type{0.25}, + d_rows.get_data()); + + GKO_ASSERT_ARRAY_EQ(d_rows, rows); +} + + +TEST_F(Pmis, ComputeStrongDepIsEquivalentToRef) +{ + initialize_data(); + auto num = system_mtx->get_size()[0]; + gko::array strong_col_idxs(ref, row_ptrs.get_const_data()[num]); + auto strong_dep = gko::matrix::SparsityCsr::create( + ref, system_mtx->get_size(), std::move(strong_col_idxs), + std::move(row_ptrs)); + gko::array d_maxabs(exec, maxabs); + auto d_strong_dep = gko::clone(exec, strong_dep); + + gko::kernels::reference::pmis::compute_strong_dep( + ref, system_mtx.get(), maxabs.get_const_data(), real_type{0.25}, + strong_dep.get()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::compute_strong_dep( + exec, d_system_mtx.get(), d_maxabs.get_const_data(), real_type{0.25}, + d_strong_dep.get()); + + GKO_ASSERT_MTX_EQ_SPARSITY(d_strong_dep, strong_dep); +} + + +TEST_F(Pmis, CountIsEquivalentToRef) +{ + initialize_data(); + auto status = gko::test::generate_random_array( + m, std::uniform_int_distribution<>(-1, 1), rand_engine, ref); + gko::array d_status(exec, status); + gko::size_type num; + gko::size_type d_num; + + gko::kernels::reference::pmis::count(ref, m, status.get_const_data(), &num); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::count( + exec, m, d_status.get_const_data(), &d_num); + + ASSERT_EQ(d_num, num); +} + + +TEST_F(Pmis, InitializeWeightAndStatusIsEquivalentToRef) +{ + initialize_data(); + auto num = system_mtx->get_size()[0]; + auto trans_strong_dep = gko::as(strong_dep->transpose()); + auto d_trans_strong_dep = gko::clone(exec, trans_strong_dep); + gko::array weight(ref, num); + gko::array status(ref, num); + gko::array d_weight(exec, num); + gko::array d_status(exec, num); + + gko::kernels::reference::pmis::initialize_weight_and_status( + ref, trans_strong_dep.get(), weight.get_data(), status.get_data()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::initialize_weight_and_status( + exec, d_trans_strong_dep.get(), d_weight.get_data(), + d_status.get_data()); + + GKO_ASSERT_ARRAY_EQ(d_status, status); + for (int i = 0; i < m; i++) { + ASSERT_EQ( + std::floor(weight.get_const_data()[i]), + std::floor(exec->copy_val_to_host(d_weight.get_const_data() + i))); + } +} + + +TEST_F(Pmis, ClassifyIsEquivalentToRef) +{ + initialize_data(); + auto num = system_mtx->get_size()[0]; + gko::array weight(ref, num); + gko::array status(ref, num); + gko::kernels::reference::pmis::initialize_weight_and_status( + ref, trans_strong_dep.get(), weight.get_data(), status.get_data()); + gko::array d_weight(exec, weight); + gko::array d_status(exec, status); + auto d_strong_dep = gko::clone(exec, strong_dep); + gko::array new_status(ref, num); + gko::array d_new_status(exec, num); + + gko::kernels::reference::pmis::classify( + ref, weight.get_data(), strong_dep.get(), status.get_const_data(), + new_status.get_data()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::classify( + exec, d_weight.get_data(), d_strong_dep.get(), + d_status.get_const_data(), d_new_status.get_data()); + + GKO_ASSERT_ARRAY_EQ(d_new_status, new_status); +} + + +TEST_F(Pmis, DirectInterpolationRowCountIsEquivalentToRef) +{ + initialize_data(); + auto num = system_mtx->get_size()[0]; + auto d_strong_dep = gko::clone(exec, strong_dep); + gko::array d_final_status(exec, final_status); + gko::array prolong_row_count(ref, num); + gko::array d_prolong_row_count(exec, num); + + gko::kernels::reference::pmis::direct_interpolation_row_count( + ref, strong_dep.get(), final_status.get_const_data(), + prolong_row_count.get_data()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::direct_interpolation_row_count( + exec, d_strong_dep.get(), d_final_status.get_const_data(), + d_prolong_row_count.get_data()); + + GKO_ASSERT_ARRAY_EQ(d_prolong_row_count, prolong_row_count); +} + + +TEST_F(Pmis, DirectInterpolationFillIsEquivalentToRef) +{ + initialize_data(); + auto num = system_mtx->get_size()[0]; + gko::array d_maxabs(exec, maxabs); + gko::array d_coarse_map(exec, coarse_map); + gko::array d_prolong_row_ptrs(exec, prolong_row_ptrs); + auto prolong_nnz = prolong_row_ptrs.get_const_data()[num]; + gko::array prolong_col_idxs(ref, prolong_nnz); + gko::array prolong_values(ref, prolong_nnz); + gko::array d_prolong_col_idxs(exec, prolong_nnz); + gko::array d_prolong_values(exec, prolong_nnz); + + gko::kernels::reference::pmis::direct_interpolation_fill( + ref, system_mtx.get(), maxabs.get_const_data(), real_type{0.25}, + coarse_map.get_const_data(), prolong_row_ptrs.get_const_data(), + prolong_col_idxs.get_data(), prolong_values.get_data()); + gko::kernels::GKO_DEVICE_NAMESPACE::pmis::direct_interpolation_fill( + exec, d_system_mtx.get(), d_maxabs.get_const_data(), real_type{0.25}, + d_coarse_map.get_const_data(), d_prolong_row_ptrs.get_const_data(), + d_prolong_col_idxs.get_data(), d_prolong_values.get_data()); + + GKO_ASSERT_ARRAY_EQ(d_prolong_col_idxs, prolong_col_idxs); + GKO_ASSERT_ARRAY_NEAR(d_prolong_values, prolong_values, + r::value); +}