Added new API call AMGX_matrix_check_symmetry (#272)

This can be run on a single process to test whether the provided matrix is symmetric (structurally or completely). Also added test to the amgx_mpi_capi example, behind the flag -cs.

Author: mattmartineau

CMakeLists.txt

@@ -173,6 +173,9 @@ endif()
 string(TOUPPER ${CMAKE_BUILD_TYPE} UPP_BUILD_NAME)
 set(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS} ${CUDA_NVCC_FLAGS_${UPP_BUILD_NAME}})
 
+# Enable device lambdas
+set(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS} --extended-lambda)
+
 # Add errors for execution space warnings and enable NVTX ranges
 set(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS} --Werror cross-execution-space-call ${NVTXRANGE_FLAG})
 

examples/amgx_mpi_capi.c

@@ -327,6 +327,25 @@ int main(int argc, char **argv)
         errAndExit("ERROR: no linear system was specified");
     }
 
+    pidx = findParamIndex(argv, argc, "-cs");
+    if(pidx != -1) 
+    {
+        int structurally_symmetric = 0;
+        int symmetric = 0;
+
+        printf("Checking matrix symmetry\n");
+        AMGX_matrix_check_symmetry(A, &structurally_symmetric, &symmetric);
+        if(symmetric) {
+            printf("A is symmetric\n");
+        }
+        else if(structurally_symmetric) {
+            printf("A is structurally_symmetric\n");
+        }
+        else {
+            printf("A is assymetric\n");
+        }
+    }
+
     //free temporary storage
     if (partition_vector != NULL) { free(partition_vector); }
 

examples/amgx_mpi_poisson7.c

@@ -129,10 +129,6 @@ int main(int argc, char **argv)
     //versions
     int major, minor;
     char *ver, *date, *time;
-    //input matrix and rhs/solution
-    int *partition_sizes = NULL;
-    int *partition_vector = NULL;
-    int partition_vector_size = 0;
     //library handles
     AMGX_Mode mode;
     AMGX_config_handle cfg;

include/amgx_c.h

@@ -608,6 +608,11 @@ AMGX_RC AMGX_API AMGX_matrix_upload_all_global_32
  const void *diag_data,
  AMGX_distribution_handle distribution);
 
+AMGX_RC AMGX_API AMGX_check_matrix_symmetry
+(AMGX_matrix_handle mtx,
+ int* structurally_symmetric,
+ int* symmetric);
+
 /*********************************************************
  * C-API deprecated
  *********************************************************/

src/amgx_c.cu

@@ -55,6 +55,7 @@
 #include <vector.h>
 #include <thrust_wrapper.h>
 
+#include "matrix_analysis.h"
 #include "amgx_types/util.h"
 #include "amgx_types/rand.h"
 #include "amgx_c_wrappers.inl"
@@ -2356,11 +2357,8 @@ extern "C" {
     {
         std::string deprecated("The AMGX_initialize_plugins API call is deprecated and can be safely removed.\n");
 
-#ifdef AMGX_WITH_MPI
         amgx_distributed_output(deprecated.c_str(), deprecated.length());
-#else
-        amgx_output(deprecated.c_str(), deprecated.length());
-#endif
+
         return AMGX_RC_OK;
     }
 
@@ -2385,11 +2383,9 @@ extern "C" {
     AMGX_RC AMGX_API AMGX_finalize_plugins()
     {
         std::string deprecated("The AMGX_finalize_plugins API call is deprecated and can be safely removed.\n");
-#ifdef AMGX_WITH_MPI
+
         amgx_distributed_output(deprecated.c_str(), deprecated.length());
-#else
-        amgx_output(deprecated.c_str(), deprecated.length());
-#endif
+
         return AMGX_RC_OK;
     }
 
@@ -5234,6 +5230,59 @@ extern "C" {
         return AMGX_RC_OK;
     }
 
+    AMGX_RC AMGX_matrix_check_symmetry(const AMGX_matrix_handle mtx, int* structurally_symmetric, int* symmetric)
+    {
+        nvtxRange nvrf(__func__);
+
+        Resources *resources = NULL;
+        AMGX_CHECK_API_ERROR(getAMGXerror(getResourcesFromMatrixHandle(mtx, &resources)), NULL)
+
+#ifdef AMGX_WITH_MPI
+            int nranks;
+            MPI_Comm_size(MPI_COMM_WORLD, &nranks);
+            if(nranks > 1) {
+                std::string err_msg("AMGX_matrix_check_symmetry cannot yet test distributed matrices, please run on 1 rank.\n");
+                amgx_distributed_output(err_msg.c_str(), err_msg.length());
+                AMGX_CHECK_API_ERROR(AMGX_ERR_BAD_PARAMETERS, resources);    //return AMGX_RC_BAD_PARAMETERS;
+            }
+#endif
+
+        AMGX_ERROR rc = AMGX_OK;
+
+        AMGX_TRIES()
+        {
+            AMGX_Mode mode = get_mode_from(mtx);
+
+            switch (mode)
+            {
+#define AMGX_CASE_LINE(CASE) case CASE: { \
+                    typedef typename TemplateMode<CASE>::Type TConfig; \
+                    typedef CWrapHandle<AMGX_matrix_handle, Matrix<TConfig>> MatrixW; \
+                    MatrixW wrapA(mtx); \
+                    Matrix<TConfig>& A = *wrapA.wrapped(); \
+                    MatrixAnalysis<TConfig> m_ana(&A); \
+                    bool verbose = false; \
+                    bool structurally_symmetric_out = false; \
+                    bool symmetric_out = false; \
+                    m_ana.checkSymmetry(structurally_symmetric_out, symmetric_out, verbose); \
+                    *structurally_symmetric = (structurally_symmetric_out ? 1 : 0); \
+                    *symmetric = (symmetric_out ? 1 : 0); \
+                    break; \
+                }
+                AMGX_FORALL_BUILDS(AMGX_CASE_LINE)
+                AMGX_FORCOMPLEX_BUILDS(AMGX_CASE_LINE)
+#undef AMGX_CASE_LINE
+
+                default:
+                    return AMGX_RC_BAD_MODE;
+            }
+        }
+
+        AMGX_CATCHES(rc)
+        AMGX_CHECK_API_ERROR(rc, resources)
+        return AMGX_RC_OK;
+    }
+
     int AMGX_Debug_get_resource_count(AMGX_resources_handle rsc)
     {
         return ((ResourceW *)rsc)->wrapped().use_count();

src/matrix_analysis.cu

@@ -134,7 +134,74 @@ void MatrixAnalysis<T_Config>::checkSymmetry(bool &structuralSymmetric, bool &sy
 
     if (TConfig::memSpace == AMGX_device)
     {
-        FatalError("Device version not implemented", AMGX_ERR_NOT_IMPLEMENTED);
+        // choose epsilon
+        double eps = 1e-12;
+
+        if (types::PODTypes<ValueTypeA>::vec_prec == AMGX_vecFloat)
+        {
+            eps = 1e-7;
+        }
+
+        const int bx = A->get_block_dimx();
+        const int by = A->get_block_dimy();
+        const int nnz = A->get_num_nz();
+
+        amgx::thrust::device_vector<bool> symm_d(1, true);
+        amgx::thrust::device_vector<bool> struct_symm_d(1, true);
+        auto* symmetric_d = amgx::thrust::raw_pointer_cast(symm_d.data());
+        auto* structurally_symmetric_d = amgx::thrust::raw_pointer_cast(struct_symm_d.data());
+        auto* col_indices = A->col_indices.raw();
+        auto* row_offsets = A->row_offsets.raw();
+        auto* values = A->values.raw();
+
+        amgx::thrust::for_each_n(thrust::device, amgx::thrust::counting_iterator<int>(0), A->get_num_rows(), [=] __device__ (int i)
+        {
+            for (int jj = row_offsets[i]; jj < row_offsets[i + 1]; jj++)
+            {
+                // check structure exists
+                int j = col_indices[jj];
+
+                // ignore diagonal
+                if (i == j) { continue; }
+
+                // loop over row j, search for column i
+                bool found_on_row = false;
+
+                for (int kk = row_offsets[j]; kk < row_offsets[j + 1]; kk++)
+                {
+                    int k = col_indices[kk];
+
+                    if (k == i)
+                    {
+                        found_on_row = true;
+
+                        // check values
+                        // check all elements
+                        const int blocksize = bx * by;
+
+                        for (int m = 0; m < bx * by; m++)
+                        {
+                            if (types::util<ValueTypeA>::abs(values[jj * blocksize + m] - values[kk * blocksize + m]) > eps)
+                            {
+                                symmetric_d[0] = false;
+                            }
+                        }
+
+                        break;
+                    }
+                }
+
+                // if we didn't find the element, non-symmetric
+                if (!found_on_row)
+                {
+                    structurally_symmetric_d[0] = false;
+                    symmetric_d[0] = false;
+                }
+            }
+        });
+
+        symmetric = symm_d[0];
+        structuralSymmetric = struct_symm_d[0];
     }
     else if (TConfig::memSpace == AMGX_host)
     {