Add KSP and classical AMG stage test coverage

ci/test.sh

@@ -20,7 +20,5 @@ mkdir -p "${BUILD_DIR}"
     cd "${BUILD_DIR}"
     cmake ..
     make -j 8 all
-    # WIP: test_launcher is allowed to fail; not all tests pass
-    set +e
-    ./tests/amgx_tests_launcher
+    ./src/amgx_tests_launcher
 )

include/classical/strength/all.h

@@ -22,6 +22,7 @@ class Strength_All : public Strength_Base<T_Config>
         {
             return true;
         }
+        bool marks_all_connections() const { return true; }
 };
 
 template<class T_Config>

include/classical/strength/strength_base.h

@@ -42,6 +42,7 @@ class Strength_BaseBase : public Strength<T_Config>
 
         __host__ __device__
         virtual bool strongly_connected(ValueType val, ValueType threshold, ValueType diagonal) = 0;
+        virtual bool marks_all_connections() const { return false; }
     protected:
         virtual void computeStrongConnectionsAndWeights_1x1(Matrix<T_Config> &A,
                 BVector &s_con,
@@ -98,4 +99,3 @@ class  Strength_Base< TemplateConfig<AMGX_device, t_vecPrec, t_matPrec, t_indPre
 };
 
 } // namespace amgx
-

src/CMakeLists.txt

@@ -15,7 +15,9 @@ target_include_directories(amgx_libs PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../inclu
 
 set(AMGX_INCLUDES ../include)
 
-set(tests_all ${tests_all} testframework.cu test_utils.cu unit_test.cu)
+file(GLOB TESTS CONFIGURE_DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/tests/*.cu")
+
+set(tests_all ${TESTS} testframework.cu test_utils.cu unit_test.cu)
 
 add_library(amgx_tests_libs OBJECT ${tests_all})
 
@@ -106,4 +108,3 @@ install(FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/configs/eigen_configs/POWER_ITERATION"
     "${CMAKE_CURRENT_SOURCE_DIR}/configs/eigen_configs/SUBSPACE_ITERATION"
     DESTINATION "lib/configs/eigen_configs")
-

src/classical/strength/strength_base.cu

@@ -189,7 +189,8 @@ void computeStrongConnectionsAndWeightsKernel( const IndexType *A_rows,
         ValueType alpha,
         ValueType *row_sum,
         const double max_row_sum,
-        int64_t base_index)
+        int64_t base_index,
+        bool mark_all_connections)
 {
     // One warp works on each row and hence one iteration handles
     // num_warps*numBlock rows. This means atomicAdd() is inevitable.
@@ -318,7 +319,8 @@ void computeStrongConnectionsAndWeightsKernel( const IndexType *A_rows,
             {
                 bool is_off_diagonal = aRowIt < aRowEnd && aColId != aRowId;
                 is_strongly_connected = is_off_diagonal &&
-                                        stronglyConnectedAHat( aValue, s_threshold[warpId], s_diag[warpId] );
+                                        (mark_all_connections ||
+                                         stronglyConnectedAHat( aValue, s_threshold[warpId], s_diag[warpId] ));
             }
 
             if ( is_strongly_connected && aRowIt < aRowEnd && aColId < A_num_rows)
@@ -348,7 +350,8 @@ void computeStrongConnectionsAndWeightsKernel_opt( const IndexType *A_rows,
         ValueType alpha,
         ValueType *row_sum,
         const double max_row_sum,
-        int64_t base_index)
+        int64_t base_index,
+        bool mark_all_connections)
 {
     // One warp works on each row and hence one iteration handles
     // num_warps*numBlock rows. This means atomicAdd() is inevitable.
@@ -430,7 +433,8 @@ void computeStrongConnectionsAndWeightsKernel_opt( const IndexType *A_rows,
             {
                 bool is_off_diagonal = aColId != aRowId;
                 is_strongly_connected = is_off_diagonal &&
-                                        stronglyConnectedAHat( aValue, threshold, diag );
+                                        (mark_all_connections ||
+                                         stronglyConnectedAHat( aValue, threshold, diag ));
             }
 
             if ( is_strongly_connected && aColId < A_num_rows)
@@ -511,6 +515,7 @@ computeStrongConnectionsAndWeights_1x1(Matrix_d &A,
     bool *s_con_ptr = s_con.raw();
     float *weights_ptr = weights.raw();
     bool compute_row_sum = (max_row_sum < 1.0);
+    const bool mark_all_connections = this->marks_all_connections();
 
     if (A.get_num_rows() == 0) { compute_row_sum = false; }
 
@@ -548,7 +553,8 @@ computeStrongConnectionsAndWeights_1x1(Matrix_d &A,
                             this->alpha,
                             compute_row_sum ? sums_ptr.raw() : NULL,
                             max_row_sum,
-                            0);
+                            0,
+                            mark_all_connections);
                 cudaCheckError();
             }
             else {
@@ -563,7 +569,8 @@ computeStrongConnectionsAndWeights_1x1(Matrix_d &A,
                             this->alpha,
                             compute_row_sum ? sums_ptr.raw() : NULL,
                             max_row_sum,
-                            A.manager->base_index());
+                            A.manager->base_index(),
+                            mark_all_connections);
                 cudaCheckError();
             }
         }
@@ -585,7 +592,8 @@ computeStrongConnectionsAndWeights_1x1(Matrix_d &A,
                             this->alpha,
                             compute_row_sum ? sums_ptr.raw() : NULL,
                             max_row_sum,
-                            0);
+                            0,
+                            mark_all_connections);
                 cudaCheckError();
             }
             else {
@@ -600,7 +608,8 @@ computeStrongConnectionsAndWeights_1x1(Matrix_d &A,
                             this->alpha,
                             compute_row_sum ? sums_ptr.raw() : NULL,
                             max_row_sum,
-                            A.manager->base_index());
+                            A.manager->base_index(),
+                            mark_all_connections);
                 cudaCheckError();
             }
         }
@@ -711,4 +720,3 @@ AMGX_FORALL_BUILDS(AMGX_CASE_LINE)
 #undef AMGX_CASE_LINE
 
 } // namespace amgx
-

src/scalers/binormalization.cu

@@ -145,7 +145,7 @@ struct std_f
 };
 
 // scaled the matrix using diag(F)*A*diag(G), f = sqrt(fabs(x)), g = sqrt(fabs(y))
-template <typename IndexType, typename MatrixType, typename VectorType>
+template <ScaleDirection direction, typename IndexType, typename MatrixType, typename VectorType>
 __global__
 void scaleMatrixDevice(int rows, IndexType *offsets, IndexType *indices, MatrixType *values,
                        VectorType *x)
@@ -158,7 +158,16 @@ void scaleMatrixDevice(int rows, IndexType *offsets, IndexType *indices, MatrixT
         {
             int j = indices[jj];
             VectorType fj = fabs(x[j]);
-            values[jj] *= sqrt(fabs(fi * fj));
+            const VectorType scale = sqrt(fabs(fi * fj));
+
+            if (direction == SCALE)
+            {
+                values[jj] *= scale;
+            }
+            else
+            {
+                values[jj] /= scale;
+            }
         }
     }
 }
@@ -375,12 +384,23 @@ void BinormalizationScaler<TemplateConfig<AMGX_device, t_vecPrec, t_matPrec, t_i
     ValueTypeB col_min = *(amgx::thrust::min_element(colnorms.begin(), colnorms.end()));
     cudaCheckError();
     printf("Original Matrix: rowmax: %e, rowmin: %e, colmax: %e, colmin: %e\n", row_max, row_min, col_max, col_min);fflush(stdout);*/
-    scaleMatrixDevice <<< 4096, 256>>>(nrows, A.row_offsets.raw(), A.col_indices.raw(), A.values.raw(), scale_vector.raw());
-    cudaCheckError();
-    ValueTypeB C_norm = sqrt(thrust_wrapper::transform_reduce<AMGX_device>(A.values.begin(), A.values.begin() + A.get_num_nz() * A.get_block_size(), square_value<ValueTypeB>(), 0., amgx::thrust::plus<ValueTypeB>()) / nrows);
-    thrust_wrapper::transform<AMGX_device>(A.values.begin(), A.values.begin() + A.get_num_nz()*A.get_block_size(), A.values.begin(), vmul_scale_const<ValueTypeB>(1. / C_norm) );
-    thrust_wrapper::transform<AMGX_device>(scale_vector.begin(), scale_vector.end(), scale_vector.begin(), vmul_scale_const<ValueTypeB>(sqrt(1. / C_norm)) );
-    cudaCheckError();
+    if (scaleOrUnscale == SCALE)
+    {
+        scaleMatrixDevice<SCALE> <<< 4096, 256>>>(nrows, A.row_offsets.raw(), A.col_indices.raw(), A.values.raw(), scale_vector.raw());
+        cudaCheckError();
+        ValueTypeB C_norm = sqrt(thrust_wrapper::transform_reduce<AMGX_device>(A.values.begin(), A.values.begin() + A.get_num_nz() * A.get_block_size(), square_value<ValueTypeB>(), 0., amgx::thrust::plus<ValueTypeB>()) / nrows);
+        thrust_wrapper::transform<AMGX_device>(A.values.begin(), A.values.begin() + A.get_num_nz()*A.get_block_size(), A.values.begin(), vmul_scale_const<ValueTypeB>(1. / C_norm) );
+        // scale_vector stores squared row/column factors and is square-rooted
+        // by scaleMatrixDevice/scaleVector.  Fold the matrix normalization into
+        // that squared factor so UNSCALE exactly reverses the operation.
+        thrust_wrapper::transform<AMGX_device>(scale_vector.begin(), scale_vector.end(), scale_vector.begin(), vmul_scale_const<ValueTypeB>(1. / C_norm) );
+        cudaCheckError();
+    }
+    else
+    {
+        scaleMatrixDevice<UNSCALE> <<< 4096, 256>>>(nrows, A.row_offsets.raw(), A.col_indices.raw(), A.values.raw(), scale_vector.raw());
+        cudaCheckError();
+    }
     /*thrust_wrapper::fill<AMGX_device>(rownorms.begin(), rownorms.end(), 0.);
       thrust_wrapper::fill<AMGX_device>(colnorms.begin(), colnorms.end(), 0.);
     getColRowNorms<<<4096,256>>>(nrows, A.row_offsets.raw(), A.col_indices.raw(), A.values.raw(), rownorms.raw(), colnorms.raw());
@@ -391,7 +411,6 @@ void BinormalizationScaler<TemplateConfig<AMGX_device, t_vecPrec, t_matPrec, t_i
     col_min = *(amgx::thrust::min_element(colnorms.begin(), colnorms.end()));
     cudaCheckError();
     printf("Scaled Matrix: rowmax: %e, rowmin: %e, colmax: %e, colmin: %e\n", row_max, row_min, col_max, col_min);fflush(stdout);*/
-    exit(0);
 }
 
 // Setup on Host
@@ -515,4 +534,3 @@ AMGX_FORALL_BUILDS(AMGX_CASE_LINE)
 
 
 } // namespace amgx
-

src/scalers/nbinormalization.cu

@@ -493,8 +493,12 @@ void NBinormalizationScaler<TemplateConfig<AMGX_device, t_vecPrec, t_matPrec, t_
             this->norm_coef = sqrt(thrust_wrapper::transform_reduce<AMGX_device>(A.values.begin(), A.values.begin() + A.get_num_nz() * A.get_block_size(), square_value<ValueTypeB>(), 0., amgx::thrust::plus<ValueTypeB>()) / A.get_num_rows());
             cudaCheckError();
             thrust_wrapper::transform<AMGX_device>(A.values.begin(), A.values.begin() + A.get_num_nz()*A.get_block_size(), A.values.begin(), vmul_scale_const<ValueTypeB>(1. / this->norm_coef) );
-            thrust_wrapper::transform<AMGX_device>(left_scale.begin(), left_scale.end(), left_scale.begin(), vmul_scale_const<ValueTypeB>(sqrt(1. / this->norm_coef)) );
-            thrust_wrapper::transform<AMGX_device>(right_scale.begin(), right_scale.end(), right_scale.begin(), vmul_scale_const<ValueTypeB>(sqrt(1. / this->norm_coef)) );
+            // left_scale/right_scale store squared factors.  Each is
+            // square-rooted by the matrix/vector scaling kernels, so folding
+            // 1/norm_coef into both squared factors gives the required
+            // combined 1/norm_coef matrix normalization.
+            thrust_wrapper::transform<AMGX_device>(left_scale.begin(), left_scale.end(), left_scale.begin(), vmul_scale_const<ValueTypeB>(1. / this->norm_coef) );
+            thrust_wrapper::transform<AMGX_device>(right_scale.begin(), right_scale.end(), right_scale.begin(), vmul_scale_const<ValueTypeB>(1. / this->norm_coef) );
             cudaCheckError();
             /*thrust_wrapper::fill<AMGX_device>(rownorms.begin(), rownorms.end(), 0.);
               thrust_wrapper::fill<AMGX_device>(colnorms.begin(), colnorms.end(), 0.);
@@ -644,4 +648,3 @@ AMGX_FORALL_BUILDS(AMGX_CASE_LINE)
 
 
 } // namespace amgx
-

src/tests/classical_coarsening_variants.cu

@@ -0,0 +1,114 @@
+// SPDX-FileCopyrightText: 2011 - 2025 NVIDIA CORPORATION. All Rights Reserved.
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include "unit_test.h"
+#include <classical/interpolators/common.h>
+#include <classical/selectors/selector.h>
+#include <classical/strength/strength_base.h>
+#include <matrix_io.h>
+
+namespace amgx
+{
+
+DECLARE_UNITTEST_BEGIN(ClassicalCoarseningVariants);
+
+typedef Vector<typename TConfig::template setVecPrec<AMGX_vecBool>::Type> BVector;
+typedef Vector<typename TConfig::template setVecPrec<AMGX_vecFloat>::Type> FVector;
+typedef Vector<typename TConfig_h::template setVecPrec<AMGX_vecBool>::Type> BVector_h;
+
+void check_selector(const char *selector_name, const Matrix_h &host_matrix, bool check_independence = true)
+{
+    AMG_Config cfg;
+    std::string parameters = std::string("selector=") + selector_name +
+                             ", strength=AHAT, strength_threshold=0.25, "
+                             "determinism_flag=1, use_opt_kernels=0";
+    UNITTEST_ASSERT_EQUAL(cfg.parseParameterString(parameters.c_str()), AMGX_OK);
+
+    Matrix<TConfig> A = host_matrix;
+    BVector strong_connections(A.get_num_nz(), false);
+    FVector weights(A.get_num_rows(), 0.0f);
+    const bool compatible_relaxation = std::string(selector_name) == "CR";
+    IVector cf_map(A.get_num_rows(), compatible_relaxation ? FINE : UNASSIGNED);
+    IVector scratch(A.get_num_rows(), 0);
+
+    Strength<TConfig> *strength = StrengthFactory<TConfig>::allocate(cfg, "default");
+    UNITTEST_ASSERT_TRUE(strength != NULL);
+    strength->computeStrongConnectionsAndWeights(A, strong_connections, weights, 1.1);
+    delete strength;
+
+    classical::Selector<TConfig> *selector =
+        classical::SelectorFactory<TConfig>::allocate(cfg, "default");
+    UNITTEST_ASSERT_TRUE(selector != NULL);
+    selector->markCoarseFinePoints(A, weights, strong_connections, cf_map, scratch);
+    delete selector;
+
+    IVector_h host_cf_map = cf_map;
+    BVector_h host_connections = strong_connections;
+    int num_coarse = 0;
+
+    if (!check_independence) { return; }
+
+    for (int row = 0; row < host_matrix.get_num_rows(); ++row)
+    {
+        const int state = host_cf_map[row];
+        this->PrintOnFail("%s left invalid state %d at row %d", selector_name, state, row);
+        UNITTEST_ASSERT_TRUE(state == COARSE || state == FINE || state == STRONG_FINE);
+        num_coarse += state == COARSE;
+    }
+
+    this->PrintOnFail("%s selected no coarse points", selector_name);
+    UNITTEST_ASSERT_TRUE(num_coarse > 0);
+    UNITTEST_ASSERT_TRUE(num_coarse < host_matrix.get_num_rows());
+
+    for (int row = 0; row < host_matrix.get_num_rows(); ++row)
+    {
+        if (host_cf_map[row] != COARSE) { continue; }
+
+        for (int jj = host_matrix.row_offsets[row]; jj < host_matrix.row_offsets[row + 1]; ++jj)
+        {
+            const int col = host_matrix.col_indices[jj];
+            if (col != row && host_connections[jj])
+            {
+                this->PrintOnFail("%s selected strongly connected coarse rows %d and %d",
+                                  selector_name, row, col);
+                UNITTEST_ASSERT_TRUE(host_cf_map[col] != COARSE);
+            }
+        }
+    }
+}
+
+void run()
+{
+    Matrix_h A;
+    A.set_initialized(0);
+    A.addProps(CSR);
+    MatrixCusp<TConfig_h, cusp::csr_format> wrapped_A(&A);
+    cusp::gallery::poisson5pt(wrapped_A, 8, 8);
+    A.computeDiagonal();
+    A.set_initialized(1);
+
+    const char *selectors[] = {"PMIS", "HMIS"};
+    for (const char *selector : selectors) { check_selector(selector, A); }
+
+    if (TConfig::memSpace == AMGX_device)
+    {
+        check_selector("AGGRESSIVE_PMIS", A);
+        check_selector("AGGRESSIVE_HMIS", A);
+        // Compatible relaxation is device-only and does not construct a
+        // maximal independent set, so validate its complete C/F assignment
+        // without imposing the PMIS/HMIS independence invariant.
+        check_selector("CR", A, false);
+    }
+}
+
+DECLARE_UNITTEST_END(ClassicalCoarseningVariants);
+
+ClassicalCoarseningVariants<TemplateMode<AMGX_mode_hDDI>::Type> ClassicalCoarseningVariants_hDDI;
+ClassicalCoarseningVariants<TemplateMode<AMGX_mode_hDFI>::Type> ClassicalCoarseningVariants_hDFI;
+ClassicalCoarseningVariants<TemplateMode<AMGX_mode_hFFI>::Type> ClassicalCoarseningVariants_hFFI;
+ClassicalCoarseningVariants<TemplateMode<AMGX_mode_dDDI>::Type> ClassicalCoarseningVariants_dDDI;
+ClassicalCoarseningVariants<TemplateMode<AMGX_mode_dDFI>::Type> ClassicalCoarseningVariants_dDFI;
+ClassicalCoarseningVariants<TemplateMode<AMGX_mode_dFFI>::Type> ClassicalCoarseningVariants_dFFI;
+
+} // namespace amgx