Add AMGX preconditioner support with GPU types and introduce backend agnostic AmgxInterface

CMakeLists.txt

@@ -433,6 +433,7 @@ macro (tests_hook)
     endif()
   endif()
   if(AMGX_FOUND)
+    list(APPEND tests_SOURCES tests/gpuistl/test_AmgxInterface.cpp)
     list(APPEND tests_SOURCES tests/test_AmgxPreconditioner.cpp)
   endif()
 endmacro (tests_hook)
@@ -790,6 +791,7 @@ if(HYPRE_FOUND)
 endif()
 
 if(AMGX_FOUND)
+  set_tests_properties(AmgxInterface PROPERTIES LABELS gpu_cuda)
   set_tests_properties(AmgxPreconditioner PROPERTIES LABELS gpu_cuda)
 endif()
 

CMakeLists_files.cmake

@@ -336,6 +336,7 @@ if (HAVE_CUDA)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/CuSparseResource.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/CuSparseResource_impl.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/safe_conversion.hpp)
+  ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/gpu_type_detection.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/cublas_wrapper.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/cusparse_wrapper.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/cusparse_constants.hpp)
@@ -358,8 +359,8 @@ if (HAVE_CUDA)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg PreconditionerCPUMatrixToGPUMatrix.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg ISTLSolverGPUISTL.hpp)
   ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/FlexibleSolverWrapper.hpp)
+  ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg AmgxInterface.hpp)
 
-  
   if(MPI_FOUND)
     ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg GpuOwnerOverlapCopy.hpp)
     ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg GpuSender.hpp)

opm/simulators/linalg/AmgxPreconditioner.hpp

@@ -25,15 +25,17 @@
 #include <opm/common/TimingMacros.hpp>
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/gpuistl/AmgxInterface.hpp>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
 
 #include <amgx_c.h>
 
-#include <vector>
+namespace Amgx
+{
 
-namespace Amgx {
+using AmgxInterface = Opm::gpuistl::AmgxInterface;
 
 /**
  * @brief Configuration structure for AMGX parameters.
@@ -54,7 +56,8 @@ struct AmgxConfig {
     double strength_threshold = 0.5;
     int max_iters = 1;
 
-    explicit AmgxConfig(const Opm::PropertyTree& prm) {
+    explicit AmgxConfig(const Opm::PropertyTree& prm)
+    {
         determinism_flag = prm.get<int>("determinism_flag", determinism_flag);
         print_grid_stats = prm.get<int>("print_grid_stats", print_grid_stats);
         print_solve_stats = prm.get<int>("print_solve_stats", print_solve_stats);
@@ -91,14 +94,16 @@ struct AmgxConfig {
  *
  * This class provides an interface to the AMG preconditioner from the AMGX library.
  * It is designed to work with matrices, update vectors, and defect vectors specified
- * by the template parameters.
+ * by the template parameters. The AmgxInterface class provides a unified interface
+ * to AMGX's functionality, allowing for easy switching between CPU and GPU input
+ * data types.
  *
- * @tparam M The matrix type the preconditioner is for.
- * @tparam X The type of the update vector.
- * @tparam Y The type of the defect vector.
+ * @tparam M The matrix type
+ * @tparam X The vector type for the solution
+ * @tparam Y The vector type for the right-hand side
  */
-template<class M, class X, class Y>
-class AmgxPreconditioner : public Dune::PreconditionerWithUpdate<X,Y>
+template <class M, class X, class Y>
+class AmgxPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
 {
 public:
     //! \brief The matrix type the preconditioner is for
@@ -115,88 +120,59 @@ public:
     static constexpr int block_size = 1;
 
     /**
-     * @brief Constructor for the AmgxPreconditioner class.
+     * @brief Constructor for AmgxPreconditioner
      *
-     * Initializes the preconditioner with the given matrix and property tree.
+     * Initializes the AMGX preconditioner with the given matrix and property tree.
      *
-     * @param A The matrix for which the preconditioner is constructed.
-     * @param prm The property tree containing configuration parameters.
+     * @param A The matrix for which to construct the preconditioner
+     * @param prm The property tree containing configuration parameters
      */
     AmgxPreconditioner(const M& A, const Opm::PropertyTree prm)
         : A_(A)
-        , N_(A.N())
-        , nnz_(A.nonzeroes())
+        , setup_frequency_(prm.get<int>("setup_frequency", 30))
+        , update_counter_(0)
     {
         OPM_TIMEBLOCK(prec_construct);
 
         // Create configuration
         AmgxConfig config(prm);
-        AMGX_SAFE_CALL(AMGX_config_create(&cfg_, config.toString().c_str()));
-        AMGX_SAFE_CALL(AMGX_resources_create_simple(&rsrc_, cfg_));
+        cfg_ = AmgxInterface::createConfig(config.toString());
+        rsrc_ = AmgxInterface::createResources(cfg_);
 
-        // Setup frequency is set in the property tree
-        setup_frequency_ = prm.get<int>("setup_frequency", 30);
+        // Determine appropriate AMGX mode based on matrix and vector types
+        amgx_mode_ = AmgxInterface::determineAmgxMode<matrix_field_type, vector_field_type>();
 
-        // Select appropriate AMGX mode based on matrix and vector scalar types
-        AMGX_Mode amgx_mode;
-        if constexpr (std::is_same_v<matrix_field_type, double> && std::is_same_v<vector_field_type, double>) {
-            amgx_mode = AMGX_mode_dDDI;
-        } else if constexpr (std::is_same_v<matrix_field_type, float> && std::is_same_v<vector_field_type, double>) {
-            amgx_mode = AMGX_mode_dDFI;
-        } else if constexpr (std::is_same_v<matrix_field_type, float> && std::is_same_v<vector_field_type, float>) {
-            amgx_mode = AMGX_mode_dFFI;
-        } else {
-            OPM_THROW(std::runtime_error, "Unsupported combination of matrix and vector types in AmgxPreconditioner");
-        }
+        // Create solver, matrix, and vector handles, given the mode
+        solver_ = AmgxInterface::createSolver(rsrc_, amgx_mode_, cfg_);
+        A_amgx_ = AmgxInterface::createMatrix(rsrc_, amgx_mode_);
+        x_amgx_ = AmgxInterface::createVector(rsrc_, amgx_mode_);
+        b_amgx_ = AmgxInterface::createVector(rsrc_, amgx_mode_);
 
-        // Create solver and matrix/vector handles with selected mode
-        AMGX_SAFE_CALL(AMGX_solver_create(&solver_, rsrc_, amgx_mode, cfg_));
-        AMGX_SAFE_CALL(AMGX_matrix_create(&A_amgx_, rsrc_, amgx_mode));
-        AMGX_SAFE_CALL(AMGX_vector_create(&x_amgx_, rsrc_, amgx_mode));
-        AMGX_SAFE_CALL(AMGX_vector_create(&b_amgx_, rsrc_, amgx_mode));
+        // Initialize matrix structure and values
+        AmgxInterface::initializeMatrix(A_, A_amgx_);
 
-        // Setup matrix structure
-        std::vector<int> row_ptrs(N_ + 1);
-        std::vector<int> col_indices(nnz_);
-        setupSparsityPattern(row_ptrs, col_indices);
+        // Initialize vectors with proper dimensions
+        const int N = A_.N();
+        AmgxInterface::initializeVector(N, block_size, x_amgx_);
+        AmgxInterface::initializeVector(N, block_size, b_amgx_);
 
-        // initialize matrix with values
-        const matrix_field_type* values = &(A_[0][0][0][0]);
-        AMGX_SAFE_CALL(AMGX_pin_memory(const_cast<matrix_field_type*>(values), sizeof(matrix_field_type) * nnz_ * block_size * block_size));
-        AMGX_SAFE_CALL(AMGX_matrix_upload_all(A_amgx_, N_, nnz_, block_size, block_size,
-                                             row_ptrs.data(), col_indices.data(),
-                                             values, nullptr));
+        // Perform initial update
         update();
     }
 
     /**
-     * @brief Destructor for the AmgxPreconditioner class.
+     * @brief Destructor for AmgxPreconditioner
      *
      * Cleans up resources allocated by the preconditioner.
      */
     ~AmgxPreconditioner()
     {
-        const matrix_field_type* values = &(A_[0][0][0][0]);
-        AMGX_SAFE_CALL(AMGX_unpin_memory(const_cast<matrix_field_type*>(values)));
-        if (solver_) {
-            AMGX_SAFE_CALL(AMGX_solver_destroy(solver_));
-        }
-        if (x_amgx_) {
-            AMGX_SAFE_CALL(AMGX_vector_destroy(x_amgx_));
-        }
-        if (b_amgx_) {
-            AMGX_SAFE_CALL(AMGX_vector_destroy(b_amgx_));
-        }
-        if (A_amgx_) {
-            AMGX_SAFE_CALL(AMGX_matrix_destroy(A_amgx_));
-        }
-        // Destroying resources and config crashes when reinitializing
-        //if (rsrc_) {
-        //    AMGX_SAFE_CALL(AMGX_resources_destroy(rsrc_));
-        //}
-        //if (cfg_) {
-        //    AMGX_SAFE_CALL(AMGX_config_destroy(cfg_));
-        //}
+        AmgxInterface::destroySolver(solver_);
+        AmgxInterface::destroyVector(x_amgx_);
+        AmgxInterface::destroyVector(b_amgx_);
+        AmgxInterface::destroyMatrix(A_amgx_, A_);
+        AmgxInterface::destroyResources(rsrc_);
+        AmgxInterface::destroyConfig(cfg_);
     }
 
     /**
@@ -207,14 +183,16 @@ public:
      * @param v The update vector.
      * @param d The defect vector.
      */
-    void pre(X& /*v*/, Y& /*d*/) override {
+    void pre(X& /*v*/, Y& /*d*/) override
+    {
     }
 
     /**
      * @brief Applies the preconditioner to a vector.
      *
      * Performs one AMG cycle to solve the system.
-     * Involves uploading vectors to AMGX, applying the preconditioner, and downloading the result.
+     * Involves uploading vectors to AMGX, applying the preconditioner,
+     * and transferring the result back to the vector.
      *
      * @param v The update vector.
      * @param d The defect vector.
@@ -223,15 +201,15 @@ public:
     {
         OPM_TIMEBLOCK(prec_apply);
 
-        // Upload vectors to AMGX
-        AMGX_SAFE_CALL(AMGX_vector_upload(x_amgx_, N_, block_size, &v[0][0]));
-        AMGX_SAFE_CALL(AMGX_vector_upload(b_amgx_, N_, block_size, &d[0][0]));
+        // Transfer vectors to AMGX
+        AmgxInterface::transferVectorToAmgx(v, x_amgx_);
+        AmgxInterface::transferVectorToAmgx(d, b_amgx_);
 
         // Apply preconditioner
         AMGX_SAFE_CALL(AMGX_solver_solve(solver_, b_amgx_, x_amgx_));
 
-        // Download result
-        AMGX_SAFE_CALL(AMGX_vector_download(x_amgx_, &v[0][0]));
+        // Transfer result back
+        AmgxInterface::transferVectorFromAmgx(x_amgx_, v);
     }
 
     /**
@@ -241,7 +219,8 @@ public:
      *
      * @param v The update vector.
      */
-    void post(X& /*v*/) override {
+    void post(X& /*v*/) override
+    {
     }
 
     /**
@@ -252,7 +231,10 @@ public:
     void update() override
     {
         OPM_TIMEBLOCK(prec_update);
-        copyMatrixToAmgx();
+
+        AmgxInterface::updateMatrixValues(A_, A_amgx_);
+
+        // Setup or resetup the solver based on update counter
         if (update_counter_ == 0) {
             AMGX_SAFE_CALL(AMGX_solver_setup(solver_, A_amgx_));
         } else {
@@ -268,7 +250,7 @@ public:
     /**
      * @brief Returns the solver category.
      *
-     * @return The solver category, which is sequential.
+     * @return The solver category (sequential).
      */
     Dune::SolverCategory::Category category() const override
     {
@@ -278,69 +260,32 @@ public:
     /**
      * @brief Checks if the preconditioner has a perfect update.
      *
-     * @return True, indicating that the preconditioner can be perfectly updated.
+     * @return True, as the preconditioner can be perfectly updated.
      */
     bool hasPerfectUpdate() const override
     {
         // The AMG hierarchy of the Amgx preconditioner can depend on the values of the matrix, so it must be recreated
-        // when the matrix values change, at given frequency. Since this is handled internally, we return true.
+        // when the matrix values change, at a given frequency. Since this is handled internally, we return true.
         return true;
     }
 
 private:
-    /**
-     * @brief Sets up the sparsity pattern for the AMGX matrix.
-     *
-     * This method initializes the row pointers and column indices for the AMGX matrix.
-     *
-     * @param row_ptrs The row pointers for the AMGX matrix.
-     * @param col_indices The column indices for the AMGX matrix.
-     */
-    void setupSparsityPattern(std::vector<int>& row_ptrs, std::vector<int>& col_indices)
-    {
-        int pos = 0;
-        row_ptrs[0] = 0;
-        for (auto row = A_.begin(); row != A_.end(); ++row) {
-            for (auto col = row->begin(); col != row->end(); ++col) {
-                col_indices[pos++] = col.index();
-            }
-            row_ptrs[row.index() + 1] = pos;
-        }
-    }
+    // Reference to the input matrix
+    const M& A_;
 
-    /**
-     * @brief Copies the matrix values to the AMGX matrix.
-     *
-     * This method updates the AMGX matrix with the current matrix values.
-     * The method assumes that the sparsity structure is the same and that 
-     * the values are stored in a contiguous array.
-     */
-    void copyMatrixToAmgx()
-    {
-        // Get direct pointer to matrix values
-        const matrix_field_type* values = &(A_[0][0][0][0]);
-        // Indexing explanation:
-        // A_[0]             - First row of the matrix
-        //     [0]           - First block in that row
-        //        [0]        - First row within the 1x1 block
-        //           [0]     - First column within the 1x1 block
-        // update matrix with new values, assuming the sparsity structure is the same
-        AMGX_SAFE_CALL(AMGX_matrix_replace_coefficients(A_amgx_, N_, nnz_, values, nullptr));
-    }
+    // AMGX handles
+    AMGX_config_handle cfg_ = nullptr;
+    AMGX_resources_handle rsrc_ = nullptr;
+    AMGX_solver_handle solver_ = nullptr;
+    AMGX_matrix_handle A_amgx_ = nullptr;
+    AMGX_vector_handle x_amgx_ = nullptr;
+    AMGX_vector_handle b_amgx_ = nullptr;
+    AMGX_Mode amgx_mode_;
 
-    const M& A_; //!< The matrix for which the preconditioner is constructed.
-    const int N_; //!< Number of rows in the matrix.
-    const int nnz_; //!< Number of non-zero elements in the matrix.
-    // Internal variables to control AMGX setup and reuse frequency
-    int setup_frequency_ = -1; //!< Frequency of updating the AMG hierarchy
-    int update_counter_ = 0; //!< Counter for setup updates.
-
-    AMGX_config_handle cfg_ = nullptr; //!< The AMGX configuration handle.
-    AMGX_resources_handle rsrc_ = nullptr; //!< The AMGX resources handle.
-    AMGX_solver_handle solver_ = nullptr; //!< The AMGX solver handle.
-    AMGX_matrix_handle A_amgx_ = nullptr; //!< The AMGX matrix handle.
-    AMGX_vector_handle x_amgx_ = nullptr; //!< The AMGX solution vector handle.
-    AMGX_vector_handle b_amgx_ = nullptr; //!< The AMGX right-hand side vector handle.
+    // Frequency of updating the AMG hierarchy
+    int setup_frequency_;
+    // Counter for setup updates
+    int update_counter_;
 };
 
 } // namespace Amgx

opm/simulators/linalg/StandardPreconditioners_gpu_serial.hpp

@@ -89,8 +89,28 @@ struct StandardPreconditioners<Operator,
                 return std::make_shared<GPUDILU>(op.getmat(), cpuMatrix, split_matrix, tune_gpu_kernels, mixed_precision_scheme, reorder);
             });
         });
+
+        // Only add AMG preconditioners to the factory if the operator
+        // is an actual matrix operator.
+        if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
+#if HAVE_AMGX
+            // Only add AMGX for scalar matrices
+            F::addCreator("amgx", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+                // Only create AMGX preconditioner for scalar matrices
+                if (op.getmat().blockSize() == 1) {
+                    auto prm_copy = prm;
+                    prm_copy.put("setup_frequency", Opm::Parameters::Get<Opm::Parameters::CprReuseInterval>());
+                    return std::make_shared<Amgx::AmgxPreconditioner<M, V, V>>(op.getmat(), prm_copy);
+                } else {
+                    OPM_THROW(std::logic_error, "AMGX preconditioner only works with scalar matrices (block size 1)");
+                }
+            });
+#endif
+
+        }
     }
 
+
 private:
 
     /**

opm/simulators/linalg/gpuistl/AmgxInterface.hpp

@@ -0,0 +1,631 @@
+/*
+  Copyright 2025 Equinor ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_AMGX_INTERFACE_HPP
+#define OPM_AMGX_INTERFACE_HPP
+
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpu_safe_call.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpu_type_detection.hpp>
+
+#include <amgx_c.h>
+#include <cuda_runtime.h>
+
+#include <fmt/core.h>
+
+#include <string>
+#include <vector>
+
+namespace Opm::gpuistl
+{
+
+/**
+ * @brief Exception class for AMGX errors
+ */
+class AmgxError : public std::runtime_error
+{
+public:
+    explicit AmgxError(const std::string& msg)
+        : std::runtime_error(msg)
+    {
+    }
+};
+
+/**
+ * @brief Get a descriptive error message for an AMGX error code
+ *
+ * @param err The AMGX error code
+ * @param expression The AMGX expression that caused the error
+ * @param file The source file where the error occurred
+ * @param function The function where the error occurred
+ * @param line The line number where the error occurred
+ * @return std::string The formatted error message
+ */
+inline std::string
+getAmgxErrorMessage(
+    AMGX_RC err, const std::string& expression, const std::string& file, const std::string& function, int line)
+{
+    char amgx_err_msg[4096];
+    AMGX_get_error_string(err, amgx_err_msg, sizeof(amgx_err_msg));
+
+    return fmt::format(
+        "AMGX error in expression: {}\nError code: {}\nError message: {}\nLocation: {}:{} in function {}",
+        expression,
+        static_cast<int>(err),
+        amgx_err_msg,
+        file,
+        line,
+        function);
+}
+
+/**
+ * @brief Safe call wrapper for AMGX functions
+ *
+ * Checks the return code from AMGX functions and throws an AmgxError if an error occurred.
+ *
+ * @param rc The AMGX return code to check
+ * @param expression The expression being evaluated (for error reporting)
+ * @param file The source file (typically __FILE__)
+ * @param function The function name (typically __func__)
+ * @param line The line number (typically __LINE__)
+ * @throws AmgxError if the return code indicates an error
+ */
+inline void
+amgxSafeCall(AMGX_RC rc, const std::string& expression, const std::string& file, const std::string& function, int line)
+{
+    if (rc != AMGX_RC_OK) {
+        throw AmgxError(getAmgxErrorMessage(rc, expression, file, function, line));
+    }
+}
+
+/**
+ * @brief Macro to wrap AMGX function calls with error checking
+ *
+ * Example usage:
+ * @code
+ * OPM_AMGX_SAFE_CALL(AMGX_initialize());
+ * @endcode
+ */
+#define OPM_AMGX_SAFE_CALL(expr) ::Opm::gpuistl::amgxSafeCall((expr), #expr, __FILE__, __func__, __LINE__)
+
+
+/**
+ * @brief Unified interface for AMGX operations with both CPU and GPU data structures
+ *
+ * This class provides utilities for working with AMGX resources and transferring
+ * data between CPU/GPU data structures and AMGX handles. It handles type detection
+ * and automatically chooses the most efficient transfer method:
+ * - For GPU types: Device-to-device transfers
+ * - For CPU types: Host-to-device transfers with memory pinning
+ *
+ * @note Error handling: All methods throw AmgxError exceptions when AMGX operations fail.
+ */
+class AmgxInterface
+{
+public:
+    /**
+     * @brief Initialize the AMGX library
+     *
+     * This should be called once at the start of the program before using any AMGX functionality.
+     * @throws AmgxError if initialization fails
+     */
+    static void initialize()
+    {
+        OPM_AMGX_SAFE_CALL(AMGX_initialize());
+    }
+
+    /**
+     * @brief Finalize the AMGX library
+     *
+     * This should be called once at the end of the program to release AMGX resources.
+     * @throws AmgxError if finalization fails
+     */
+    static void finalize()
+    {
+        OPM_AMGX_SAFE_CALL(AMGX_finalize());
+    }
+
+    /**
+     * @brief Create an AMGX config handle from a configuration string
+     *
+     * @param config_string Configuration string for AMGX
+     * @return AMGX_config_handle The created config handle
+     * @throws AmgxError if config creation fails
+     */
+    static AMGX_config_handle createConfig(const std::string& config_string)
+    {
+        AMGX_config_handle config;
+        OPM_AMGX_SAFE_CALL(AMGX_config_create(&config, config_string.c_str()));
+        return config;
+    }
+
+    /**
+     * @brief Create AMGX resources from a config
+     *
+     * @param config The AMGX config handle
+     * @return AMGX_resources_handle The created resources handle
+     * @throws AmgxError if resource creation fails
+     */
+    static AMGX_resources_handle createResources(AMGX_config_handle config)
+    {
+        AMGX_resources_handle resources;
+        OPM_AMGX_SAFE_CALL(AMGX_resources_create_simple(&resources, config));
+        return resources;
+    }
+
+    /**
+     * @brief Create an AMGX solver
+     *
+     * @param resources AMGX resources handle
+     * @param mode AMGX mode (precision configuration)
+     * @param config AMGX config handle
+     * @return AMGX_solver_handle The created solver handle
+     * @throws AmgxError if solver creation fails
+     */
+    static AMGX_solver_handle createSolver(AMGX_resources_handle resources, AMGX_Mode mode, AMGX_config_handle config)
+    {
+        AMGX_solver_handle solver;
+        OPM_AMGX_SAFE_CALL(AMGX_solver_create(&solver, resources, mode, config));
+        return solver;
+    }
+
+    /**
+     * @brief Create an AMGX matrix
+     *
+     * @param resources AMGX resources handle
+     * @param mode AMGX mode (precision configuration)
+     * @return AMGX_matrix_handle The created matrix handle
+     * @throws AmgxError if matrix creation fails
+     */
+    static AMGX_matrix_handle createMatrix(AMGX_resources_handle resources, AMGX_Mode mode)
+    {
+        AMGX_matrix_handle matrix;
+        OPM_AMGX_SAFE_CALL(AMGX_matrix_create(&matrix, resources, mode));
+        return matrix;
+    }
+
+    /**
+     * @brief Create an AMGX vector
+     *
+     * @param resources AMGX resources handle
+     * @param mode AMGX mode (precision configuration)
+     * @return AMGX_vector_handle The created vector handle
+     * @throws AmgxError if vector creation fails
+     */
+    static AMGX_vector_handle createVector(AMGX_resources_handle resources, AMGX_Mode mode)
+    {
+        AMGX_vector_handle vector;
+        OPM_AMGX_SAFE_CALL(AMGX_vector_create(&vector, resources, mode));
+        return vector;
+    }
+
+    /**
+     * @brief Destroy an AMGX config handle
+     *
+     * @param config The config handle to destroy
+     * @throws AmgxError if config destruction fails
+     */
+    static void destroyConfig(AMGX_config_handle config)
+    {
+        if (config) {
+            OPM_AMGX_SAFE_CALL(AMGX_config_destroy(config));
+        }
+    }
+
+    /**
+     * @brief Destroy an AMGX resources handle
+     *
+     * @param resources The resources handle to destroy
+     * @throws AmgxError if resource destruction fails
+     *
+     * @note There is a known issue where destroying resources and then
+     * reinitializing AMGX can cause crashes in some versions of the library.
+     * If you encounter this issue, consider keeping the resources alive for
+     * the entire program lifetime.
+     */
+    static void destroyResources(AMGX_resources_handle resources)
+    {
+        if (resources) {
+            OPM_AMGX_SAFE_CALL(AMGX_resources_destroy(resources));
+        }
+    }
+
+    /**
+     * @brief Destroy an AMGX solver handle
+     *
+     * @param solver The solver handle to destroy
+     * @throws AmgxError if solver destruction fails
+     */
+    static void destroySolver(AMGX_solver_handle solver)
+    {
+        if (solver) {
+            OPM_AMGX_SAFE_CALL(AMGX_solver_destroy(solver));
+        }
+    }
+
+    /**
+     * @brief Destroy an AMGX matrix handle
+     *
+     * @param matrix The matrix handle to destroy
+     * @throws AmgxError if matrix destruction fails
+     */
+    template <typename MatrixType>
+    static void destroyMatrix(AMGX_matrix_handle amgx_matrix, const MatrixType& matrix)
+    {
+        if constexpr (!is_gpu_type<MatrixType>::value) {
+            // Unpin memory for CPU matrices
+            using T = typename MatrixType::field_type;
+            const T* values = &(matrix[0][0][0][0]);
+            OPM_AMGX_SAFE_CALL(AMGX_unpin_memory(const_cast<T*>(values)));
+        }
+
+        if (amgx_matrix) {
+            OPM_AMGX_SAFE_CALL(AMGX_matrix_destroy(amgx_matrix));
+        }
+    }
+
+    /**
+     * @brief Destroy an AMGX vector handle
+     *
+     * @param vector The vector handle to destroy
+     * @throws AmgxError if vector destruction fails
+     */
+    static void destroyVector(AMGX_vector_handle vector)
+    {
+        if (vector) {
+            OPM_AMGX_SAFE_CALL(AMGX_vector_destroy(vector));
+        }
+    }
+
+    /**
+     * @brief Update an AMGX vector from a GpuVector (device-to-device transfer)
+     *
+     * Updates the AMGX vector with the contents of the GpuVector using direct device memory access.
+     * The AMGX vector must already be created with appropriate resources and mode.
+     *
+     * @param gpu_vec The source GpuVector
+     * @param amgx_vec The AMGX vector to update
+     * @throws AmgxError if the transfer fails or sizes mismatch
+     */
+    template <typename T>
+    static void updateAmgxFromGpuVector(const GpuVector<T>& gpu_vec, AMGX_vector_handle amgx_vec)
+    {
+        // Get vector size from AMGX to verify compatibility
+        int n, block_dim;
+        OPM_AMGX_SAFE_CALL(AMGX_vector_get_size(amgx_vec, &n, &block_dim));
+
+        if (n > 0 && static_cast<size_t>(n * block_dim) != gpu_vec.dim()) {
+            throw AmgxError(fmt::format("Vector size mismatch in updateAmgxFromGpuVector: "
+                                        "AMGX vector size {} vs. GpuVector size {}",
+                                        n * block_dim,
+                                        gpu_vec.dim()));
+        }
+
+        // Get raw device pointer directly from GpuVector
+        const T* device_ptr = gpu_vec.data();
+
+        // Update AMGX vector directly with the device pointer
+        OPM_AMGX_SAFE_CALL(AMGX_vector_upload(amgx_vec, n, block_dim, const_cast<T*>(device_ptr)));
+    }
+
+    /**
+     * @brief Update a GpuVector from an AMGX vector (device-to-device transfer)
+     *
+     * Updates the GpuVector with the contents of the AMGX vector using direct device memory access.
+     *
+     * @param amgx_vec The source AMGX vector
+     * @param gpu_vec The GpuVector to update
+     * @throws AmgxError if the transfer fails or sizes mismatch
+     */
+    template <typename T>
+    static void updateGpuVectorFromAmgx(AMGX_vector_handle amgx_vec, GpuVector<T>& gpu_vec)
+    {
+        // Get vector size from AMGX to verify compatibility
+        int n, block_dim;
+        OPM_AMGX_SAFE_CALL(AMGX_vector_get_size(amgx_vec, &n, &block_dim));
+
+        if (static_cast<size_t>(n * block_dim) != gpu_vec.dim()) {
+            throw AmgxError(fmt::format("Vector size mismatch in updateGpuVectorFromAmgx: "
+                                        "AMGX vector size {} vs. GpuVector size {}",
+                                        n * block_dim,
+                                        gpu_vec.dim()));
+        }
+
+        // Get destination device pointer from GpuVector
+        T* dst_device_ptr = gpu_vec.data();
+
+        // Download data directly from AMGX vector to GpuVector's device memory
+        OPM_AMGX_SAFE_CALL(AMGX_vector_download(amgx_vec, dst_device_ptr));
+    }
+
+    /**
+     * @brief Transfer vector to AMGX from any vector type (CPU or GPU)
+     *
+     * Selects the transfer method based on the vector type.
+     *
+     * @param vec Source vector (typically BlockVector (CPU) or GpuVector (GPU))
+     * @param amgx_vec Destination AMGX vector
+     * @throws AmgxError if the transfer fails
+     */
+    template <typename VectorType>
+    static void transferVectorToAmgx(const VectorType& vec, AMGX_vector_handle amgx_vec)
+    {
+        if constexpr (is_gpu_type<VectorType>::value) {
+            // GPU implementation - device-to-device transfer
+            updateAmgxFromGpuVector(vec, amgx_vec);
+        } else {
+            // CPU implementation - host-to-device transfer
+            // For BlockVector, get dimensions
+            const int N = vec.size();
+            const int block_size = 1; // Assumes scalar vectors for CPU case
+
+            // Upload vector to AMGX
+            OPM_AMGX_SAFE_CALL(AMGX_vector_upload(amgx_vec, N, block_size, &vec[0][0]));
+        }
+    }
+
+    /**
+     * @brief Transfer vector from AMGX to any vector type (CPU or GPU)
+     *
+     * Selects the transfer method based on the vector type.
+     *
+     * @param amgx_vec Source AMGX vector
+     * @param vec Destination vector (typically BlockVector (CPU) or GpuVector (GPU))
+     * @throws AmgxError if the transfer fails
+     */
+    template <typename VectorType>
+    static void transferVectorFromAmgx(AMGX_vector_handle amgx_vec, VectorType& vec)
+    {
+        if constexpr (is_gpu_type<VectorType>::value) {
+            // GPU implementation - device-to-device transfer
+            updateGpuVectorFromAmgx(amgx_vec, vec);
+        } else {
+            // CPU implementation - device-to-host transfer
+            OPM_AMGX_SAFE_CALL(AMGX_vector_download(amgx_vec, &vec[0][0]));
+        }
+    }
+
+    /**
+     * @brief Update an AMGX matrix from a GpuSparseMatrix (device-to-device transfer)
+     *
+     * Uploads the entire matrix structure and values from GpuSparseMatrix to AMGX.
+     *
+     * @param gpuSparseMatrix The source GpuSparseMatrix
+     * @param amgxMatrix The AMGX matrix to update
+     * @throws AmgxError if the transfer fails
+     */
+    template <typename T>
+    static void updateAmgxMatrixFromGpuSparseMatrix(const GpuSparseMatrix<T>& gpuSparseMatrix,
+                                                    AMGX_matrix_handle amgxMatrix)
+    {
+        // Get matrix dimensions and sparsity information
+        auto n = detail::to_int(gpuSparseMatrix.N());
+        auto nnz = detail::to_int(gpuSparseMatrix.nonzeroes());
+        auto block_size = detail::to_int(gpuSparseMatrix.blockSize());
+
+        // Get device pointers directly from GpuSparseMatrix
+        const T* values = gpuSparseMatrix.getNonZeroValues().data();
+        const int* row_ptrs = gpuSparseMatrix.getRowIndices().data();
+        const int* col_indices = gpuSparseMatrix.getColumnIndices().data();
+
+        // Update AMGX matrix with the device pointers
+        OPM_AMGX_SAFE_CALL(
+            AMGX_matrix_upload_all(amgxMatrix, n, nnz, block_size, block_size, row_ptrs, col_indices, values, nullptr));
+    }
+
+    /**
+     * @brief Update only the coefficient values of an AMGX matrix from a GpuSparseMatrix
+     *
+     * Updates just the coefficient values of an AMGX matrix without changing its sparsity pattern.
+     * This is more efficient when the matrix structure remains the same.
+     *
+     * @param gpuSparseMatrix The source GpuSparseMatrix with updated values
+     * @param amgxMatrix The AMGX matrix to update
+     * @throws AmgxError if the update fails
+     */
+    template <typename T>
+    static void updateAmgxMatrixCoefficientsFromGpuSparseMatrix(const GpuSparseMatrix<T>& gpuSparseMatrix,
+                                                                AMGX_matrix_handle amgxMatrix)
+    {
+        // Get matrix dimensions and sparsity information
+        auto n = detail::to_int(gpuSparseMatrix.N());
+        auto nnz = detail::to_int(gpuSparseMatrix.nonzeroes());
+
+        // Get device pointer to values directly from GpuSparseMatrix
+        const T* values = gpuSparseMatrix.getNonZeroValues().data();
+
+        // Update coefficients in AMGX matrix with the device pointer
+        OPM_AMGX_SAFE_CALL(AMGX_matrix_replace_coefficients(amgxMatrix, n, nnz, values, nullptr));
+    }
+
+    /**
+     * @brief Update a GpuSparseMatrix from an AMGX matrix (device-to-device transfer)
+     *
+     * Downloads the matrix data from AMGX and updates the GpuSparseMatrix with the values.
+     * The sparsity pattern is assumed to be identical, so only values are updated.
+     *
+     * @note This function is intended primarily for testing purposes and should not be used in performance-critical code paths,
+     *       as it may involve unnecessary device-to-device transfers and temporary allocations.
+     *
+     * @param amgxMatrix The source AMGX matrix
+     * @param gpuSparseMatrix The GpuSparseMatrix to update
+     * @throws AmgxError if the transfer fails
+     */
+    template <typename T>
+    static void updateGpuSparseMatrixFromAmgxMatrix(AMGX_matrix_handle amgxMatrix, GpuSparseMatrix<T>& gpuSparseMatrix)
+    {
+        // Get matrix dimensions from AMGX
+        int n, nnz, block_sizex, block_sizey;
+        OPM_AMGX_SAFE_CALL(AMGX_matrix_get_size(amgxMatrix, &n, &block_sizex, &block_sizey));
+        OPM_AMGX_SAFE_CALL(AMGX_matrix_get_nnz(amgxMatrix, &nnz));
+
+        // Allocate temporary device memory for row pointers and column indices (required by AMGX API)
+        int* temp_row_ptrs;
+        int* temp_col_indices;
+        OPM_GPU_SAFE_CALL(cudaMalloc(&temp_row_ptrs, (n + 1) * sizeof(int)));
+        OPM_GPU_SAFE_CALL(cudaMalloc(&temp_col_indices, nnz * sizeof(int)));
+
+        // Get device pointer to values in the GpuSparseMatrix
+        T* gpu_values = gpuSparseMatrix.getNonZeroValues().data();
+
+        // AMGX requires a valid pointer for diagonal data (even if unused)
+        void* diag_data_ptr = nullptr;
+
+        // Download matrix values directly from AMGX to GpuSparseMatrix
+        OPM_AMGX_SAFE_CALL(AMGX_matrix_download_all(amgxMatrix,
+                                                    temp_row_ptrs,
+                                                    temp_col_indices,
+                                                    gpu_values,
+                                                    &diag_data_ptr));
+
+        // Clean up temporary device memory
+        OPM_GPU_SAFE_CALL(cudaFree(temp_row_ptrs));
+        OPM_GPU_SAFE_CALL(cudaFree(temp_col_indices));
+    }
+
+    /**
+     * @brief Initialize an AMGX matrix from any matrix type (CPU or GPU)
+     *
+     * Selects the transfer method based on the matrix type.
+     *
+     * @param matrix Source matrix (typically BCRSMatrix (CPU) or GpuSparseMatrix (GPU))
+     * @param amgx_matrix Destination AMGX matrix
+     * @throws AmgxError if initialization fails
+     */
+    template <typename MatrixType>
+    static void initializeMatrix(const MatrixType& matrix, AMGX_matrix_handle amgx_matrix)
+    {
+        if constexpr (is_gpu_type<MatrixType>::value) {
+            // GPU implementation - device-to-device transfer
+            updateAmgxMatrixFromGpuSparseMatrix(matrix, amgx_matrix);
+        } else {
+            // CPU implementation - host-to-device transfer
+            auto N = detail::to_int(matrix.N());
+            auto nnz = detail::to_int(matrix.nonzeroes());
+            const int block_size = 1; // Assumes scalar matrices for CPU case
+
+            // Setup sparsity pattern
+            std::vector<int> row_ptrs(N + 1);
+            std::vector<int> col_indices(nnz);
+
+            int pos = 0;
+            row_ptrs[0] = 0;
+            for (auto row = matrix.begin(); row != matrix.end(); ++row) {
+                for (auto col = row->begin(); col != row->end(); ++col) {
+                    col_indices[pos++] = col.index();
+                }
+                row_ptrs[row.index() + 1] = pos;
+            }
+
+            // Get values pointer
+            using T = typename MatrixType::field_type;
+            const T* values = &(matrix[0][0][0][0]);
+            // Indexing explanation:
+            // matrix[0]             - First row of the matrix
+            //          [0]          - First block in that row
+            //            [0]        - First row within the 1x1 block
+            //              [0]      - First column within the 1x1 block
+
+            // Pin memory for CPU-to-GPU transfers
+            OPM_AMGX_SAFE_CALL(AMGX_pin_memory(const_cast<T*>(values), sizeof(T) * nnz * block_size * block_size));
+
+            // Upload to AMGX
+            OPM_AMGX_SAFE_CALL(AMGX_matrix_upload_all(amgx_matrix,
+                                                      N,
+                                                      nnz,
+                                                      block_size,
+                                                      block_size,
+                                                      row_ptrs.data(),
+                                                      col_indices.data(),
+                                                      const_cast<T*>(values),
+                                                      nullptr));
+        }
+    }
+
+    /**
+     * @brief Initialize an AMGX vector with zeros
+     *
+     * Creates a zero-initialized vector of the appropriate size.
+     *
+     * @param N Size of the vector
+     * @param block_size Block size (typically 1 for scalar vectors)
+     * @param amgx_vector The AMGX vector to initialize
+     * @throws AmgxError if initialization fails
+     */
+    static void initializeVector(int N, int block_size, AMGX_vector_handle amgx_vector)
+    {
+        OPM_AMGX_SAFE_CALL(AMGX_vector_set_zero(amgx_vector, N, block_size));
+    }
+
+    /**
+     * @brief Update matrix values in AMGX
+     *
+     * Selects the update method based on the matrix type.
+     * Updates the coefficients of an AMGX matrix, without changing its sparsity pattern.
+     *
+     * @param matrix Source matrix with updated values
+     * @param amgx_matrix AMGX matrix to update
+     * @throws AmgxError if the update fails
+     */
+    template <typename MatrixType>
+    static void updateMatrixValues(const MatrixType& matrix, AMGX_matrix_handle amgx_matrix)
+    {
+        if constexpr (is_gpu_type<MatrixType>::value) {
+            // GPU implementation - device-to-device transfer
+            updateAmgxMatrixCoefficientsFromGpuSparseMatrix(matrix, amgx_matrix);
+        } else {
+            // CPU implementation - host-to-device transfer
+            using T = typename MatrixType::field_type;
+            const T* values = &(matrix[0][0][0][0]);
+            OPM_AMGX_SAFE_CALL(AMGX_matrix_replace_coefficients(
+                amgx_matrix, matrix.N(), matrix.nonzeroes(), const_cast<T*>(values), nullptr));
+        }
+    }
+
+    /**
+     * @brief Determine the appropriate AMGX mode based on matrix and vector field types
+     *
+     * @tparam MatrixFieldType The field type of the matrix
+     * @tparam VectorFieldType The field type of the vector
+     * @return AMGX_Mode The appropriate AMGX mode
+     * @throws std::runtime_error if the type combination is not supported
+     */
+    template <typename MatrixFieldType, typename VectorFieldType>
+    static AMGX_Mode determineAmgxMode()
+    {
+        if constexpr (std::is_same_v<MatrixFieldType, double> && std::is_same_v<VectorFieldType, double>) {
+            return AMGX_mode_dDDI;
+        } else if constexpr (std::is_same_v<MatrixFieldType, float> && std::is_same_v<VectorFieldType, double>) {
+            return AMGX_mode_dDFI;
+        } else if constexpr (std::is_same_v<MatrixFieldType, float> && std::is_same_v<VectorFieldType, float>) {
+            return AMGX_mode_dFFI;
+        } else {
+            OPM_THROW(std::runtime_error,
+                      "Unsupported combination of matrix and vector types for AMGX: "
+                          + std::string(typeid(MatrixFieldType).name()) + " and "
+                          + std::string(typeid(VectorFieldType).name()));
+        }
+    }
+};
+
+} // namespace Opm::gpuistl
+
+#endif // OPM_AMGX_INTERFACE_HPP

opm/simulators/linalg/gpuistl/detail/gpu_type_detection.hpp

@@ -0,0 +1,64 @@
+/*
+  Copyright 2025 Equinor ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_GPUISTL_GPU_TYPE_DETECTION_HPP
+#define OPM_GPUISTL_GPU_TYPE_DETECTION_HPP
+
+#include <type_traits>
+
+namespace Opm::gpuistl
+{
+
+// Forward declarations
+template <typename T>
+class GpuVector;
+template <typename T>
+class GpuSparseMatrix;
+template <typename T>
+class GpuSparseMatrixGeneric;
+
+/**
+ * @brief Type trait to detect if a type is a GPU type
+ */
+template <typename T>
+struct is_gpu_type : std::false_type {
+};
+
+// Specializations for known GPU types
+template <typename T>
+struct is_gpu_type<GpuVector<T>> : std::true_type {
+};
+
+template <typename T>
+struct is_gpu_type<GpuSparseMatrix<T>> : std::true_type {
+};
+
+template <typename T>
+struct is_gpu_type<GpuSparseMatrixGeneric<T>> : std::true_type {
+};
+
+/**
+ * @brief Helper variable template for easier usage
+ */
+template <typename T>
+inline constexpr bool is_gpu_type_v = is_gpu_type<T>::value;
+
+} // namespace Opm::gpuistl
+
+#endif // OPM_GPUISTL_GPU_TYPE_DETECTION_HPP

tests/AmgxPreconditionerTestHelper.hpp

@@ -31,13 +31,18 @@
 
 #include <opm/simulators/linalg/AmgxPreconditioner.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
+#include <opm/simulators/linalg/gpuistl/detail/gpu_type_detection.hpp>
 
-template<class Matrix>
-void setupLaplace2d(int N, Matrix& mat)
+
+template <class Matrix>
+void
+setupLaplace2d(int N, Matrix& mat)
 {
-    const int nonZeroes = N*N * 5; // max 5 entries per row (diagonal + 4 neighbors)
+    const int nonZeroes = N * N * 5; // max 5 entries per row (diagonal + 4 neighbors)
     mat.setBuildMode(Matrix::row_wise);
-    mat.setSize(N*N, N*N, nonZeroes);
+    mat.setSize(N * N, N * N, nonZeroes);
 
     // Set up sparsity pattern
     for (auto row = mat.createbegin(); row != mat.createend(); ++row) {
@@ -45,15 +50,15 @@ void setupLaplace2d(int N, Matrix& mat)
         int x = i % N;
         int y = i / N;
 
-        row.insert(i);  // diagonal
-        if (x > 0)      // left neighbor
-            row.insert(i-1);
-        if (x < N-1)    // right neighbor
-            row.insert(i+1);
-        if (y > 0)      // upper neighbor
-            row.insert(i-N);
-        if (y < N-1)    // lower neighbor
-            row.insert(i+N);
+        row.insert(i); // diagonal
+        if (x > 0) // left neighbor
+            row.insert(i - 1);
+        if (x < N - 1) // right neighbor
+            row.insert(i + 1);
+        if (y > 0) // upper neighbor
+            row.insert(i - N);
+        if (y < N - 1) // lower neighbor
+            row.insert(i + N);
     }
 
     // Fill the matrix with values
@@ -66,48 +71,63 @@ void setupLaplace2d(int N, Matrix& mat)
         (*row)[i] = 4.0;
 
         // Set off-diagonal entries
-        if (x > 0)         // left neighbor
-            (*row)[i-1] = -1.0;
-        if (x < N-1)       // right neighbor
-            (*row)[i+1] = -1.0;
-        if (y > 0)         // upper neighbor
-            (*row)[i-N] = -1.0;
-        if (y < N-1)       // lower neighbor
-            (*row)[i+N] = -1.0;
+        if (x > 0) // left neighbor
+            (*row)[i - 1] = -1.0;
+        if (x < N - 1) // right neighbor
+            (*row)[i + 1] = -1.0;
+        if (y > 0) // upper neighbor
+            (*row)[i - N] = -1.0;
+        if (y < N - 1) // lower neighbor
+            (*row)[i + N] = -1.0;
     }
 }
 
-template<typename MatrixScalar, typename VectorScalar>
-void testAmgxPreconditioner()
+template <typename MatrixType>
+MatrixType
+createMatrix(int N)
+{
+    if constexpr (!Opm::gpuistl::is_gpu_type_v<MatrixType>) {
+        MatrixType matrix;
+        setupLaplace2d(N, matrix);
+        return matrix;
+    } else {
+        // For GPU matrix, create CPU matrix first then convert
+        using ScalarType = typename MatrixType::field_type;
+        using CpuMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<ScalarType, 1, 1>>;
+        CpuMatrix cpu_matrix;
+        setupLaplace2d(N, cpu_matrix);
+        return MatrixType::fromMatrix(cpu_matrix);
+    }
+}
+
+
+template <typename MatrixType, typename VectorType>
+void
+testAmgxPreconditionerGeneric()
 {
     constexpr int N = 100; // 100x100 grid
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<MatrixScalar, 1, 1>>;
-    using Vector = Dune::BlockVector<Dune::FieldVector<VectorScalar, 1>>;
 
     // Create matrix
-    Matrix matrix;
-    setupLaplace2d(N, matrix);
+    auto matrix = createMatrix<MatrixType>(N);
 
     // Create vectors
-    Vector x(N*N), b(N*N);
-    x = 100.0;    // Initial guess
-    b = 0.0;      // RHS
+    VectorType x(N * N), b(N * N);
+    x = 0.0;
+    b = 1.0;
 
     // Create operator
-    using Operator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
-    Operator op(matrix);
-
-    // Set up AMGX parameters
-    Opm::PropertyTree prm;
+    using Operator = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
+    auto op = std::make_shared<Operator>(matrix);
 
     // Create preconditioner
-    auto prec = std::make_shared<Amgx::AmgxPreconditioner<Matrix, Vector, Vector>>(matrix, prm);
+    Opm::PropertyTree prm;
+    auto prec = std::make_shared<Amgx::AmgxPreconditioner<MatrixType, VectorType, VectorType>>(matrix, prm);
 
     // Create solver
     double reduction = 1e-8;
     int maxit = 300;
     int verbosity = 0;
-    Dune::LoopSolver<Vector> solver(op, *prec, reduction, maxit, verbosity);
+    Dune::LoopSolver<VectorType> solver(*op, *prec, reduction, maxit, verbosity);
 
     // Solve
     Dune::InverseOperatorResult res;
@@ -115,7 +135,104 @@ void testAmgxPreconditioner()
 
     // Check convergence
     BOOST_CHECK(res.converged);
-    BOOST_CHECK_LT(res.reduction, 1e-8);
+    BOOST_CHECK_LT(res.reduction, reduction);
+
+    std::vector<typename VectorType::field_type> x_values;
+    if constexpr (Opm::gpuistl::is_gpu_type_v<VectorType>) {
+        x_values = x.asStdVector();
+    } else {
+        x_values.assign(&x[0][0], &x[0][0] + x.dim());
+    }
+
+    // Validate solution - check that solution has expected spatial distribution
+    double center_value = x_values[N / 2 + (N / 2) * N];
+    double near_center = x_values[(N / 2 + 1) + (N / 2) * N];
+    double far_from_center = x_values[0]; // Corner
+
+    BOOST_CHECK_GT(center_value, 0.0);
+    BOOST_CHECK_LT(near_center, center_value);
+    BOOST_CHECK_LT(far_from_center, near_center);
+}
+
+template <typename MatrixType, typename VectorType>
+void
+testAmgxPreconditionerUpdateGeneric()
+{
+    constexpr int N = 20;
+
+    // Create matrix
+    auto matrix = createMatrix<MatrixType>(N);
+
+    // Create vectors
+    VectorType x(N * N), b(N * N);
+    x = 0.0;
+    b = 1.0;
+
+    // Create preconditioner and apply it with initial values
+    Opm::PropertyTree prm;
+    auto prec = std::make_shared<Amgx::AmgxPreconditioner<MatrixType, VectorType, VectorType>>(matrix, prm);
+    prec->apply(x, b);
+
+    // Update matrix with new values and update the preconditioner
+    if constexpr (!Opm::gpuistl::is_gpu_type_v<MatrixType>) {
+        matrix *= 2.0;
+    } else {
+        matrix.getNonZeroValues() *= 2.0;
+    }
+    prec->update();
+
+    // Apply again after update but with same initial guess
+    VectorType x2(N * N);
+    x2 = 0.0;
+    prec->apply(x2, b);
+
+    std::vector<typename VectorType::field_type> x_values, x2_values;
+    if constexpr (Opm::gpuistl::is_gpu_type_v<VectorType>) {
+        x_values = x.asStdVector();
+        x2_values = x2.asStdVector();
+    } else {
+        x_values.assign(&x[0][0], &x[0][0] + x.dim());
+        x2_values.assign(&x2[0][0], &x2[0][0] + x2.dim());
+    }
+    // The vectors should be different after applying preconditioner with updated matrix
+    bool all_vectors_different = true;
+    for (size_t i = 0; i < x_values.size(); ++i) {
+        if (std::abs(x_values[i] - x2_values[i]) <= 1e-10) {
+            all_vectors_different = false;
+            break;
+        }
+    }
+    BOOST_CHECK(all_vectors_different);
+}
+
+template <typename Scalar, bool GpuInput = false>
+void
+testAmgxPreconditioner()
+{
+    if constexpr (GpuInput) {
+        using Matrix = Opm::gpuistl::GpuSparseMatrix<Scalar>;
+        using Vector = Opm::gpuistl::GpuVector<Scalar>;
+        testAmgxPreconditionerGeneric<Matrix, Vector>();
+    } else {
+        using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
+        using Vector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
+        testAmgxPreconditionerGeneric<Matrix, Vector>();
+    }
+}
+
+template <typename Scalar, bool GpuInput = false>
+void
+testAmgxPreconditionerUpdate()
+{
+    if constexpr (GpuInput) {
+        using Matrix = Opm::gpuistl::GpuSparseMatrix<Scalar>;
+        using Vector = Opm::gpuistl::GpuVector<Scalar>;
+        testAmgxPreconditionerUpdateGeneric<Matrix, Vector>();
+    } else {
+        using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
+        using Vector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
+        testAmgxPreconditionerUpdateGeneric<Matrix, Vector>();
+    }
 }
 
 #endif // TEST_AMGXPRECONDITIONER_HELPER_HPP

tests/gpuistl/test_AmgxInterface.cpp

@@ -0,0 +1,252 @@
+/*
+  Copyright 2025 Equinor ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE TestAmgxInterface
+#define BOOST_TEST_NO_MAIN
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/linalg/gpuistl/AmgxInterface.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
+#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <cstddef>
+#include <random>
+#include <vector>
+
+using namespace Opm::gpuistl;
+
+// RAII class for managing AMGX resources during tests
+class AmgxResourceGuard
+{
+public:
+    AmgxResourceGuard()
+    {
+        AmgxInterface::initialize();
+        cfg_ = AmgxInterface::createConfig("");
+        rsrc_ = AmgxInterface::createResources(cfg_);
+    }
+
+    ~AmgxResourceGuard()
+    {
+        AmgxInterface::destroyResources(rsrc_);
+        AmgxInterface::destroyConfig(cfg_);
+        AmgxInterface::finalize();
+    }
+
+    AMGX_config_handle getConfig() const
+    {
+        return cfg_;
+    }
+    AMGX_resources_handle getResources() const
+    {
+        return rsrc_;
+    }
+
+private:
+    AMGX_config_handle cfg_;
+    AMGX_resources_handle rsrc_;
+};
+
+// Test vector transfer functionality
+BOOST_AUTO_TEST_CASE(TestVectorTransfer)
+{
+    // Setup AMGX resources
+    AmgxResourceGuard resources;
+
+    // Create test data
+    const int n = 100;
+    std::vector<double> host_data(n);
+
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<double> dist(-100.0, 100.0);
+
+    // Initialize with random data
+    for (size_t i = 0; i < n; ++i) {
+        host_data[i] = dist(gen);
+    }
+
+    // Create GPU vector with test data
+    GpuVector<double> gpu_vec(host_data);
+
+    // Create AMGX vector
+    AMGX_vector_handle amgx_vec = AmgxInterface::createVector(resources.getResources(), AMGX_mode_dDDI);
+
+    // Set dimensions for the AMGX vector (needed before first transfer)
+    OPM_AMGX_SAFE_CALL(AMGX_vector_set_zero(amgx_vec, n, 1));
+
+    // Test transfer GPU -> AMGX
+    AmgxInterface::updateAmgxFromGpuVector(gpu_vec, amgx_vec);
+
+    // Create new GPU vector to receive AMGX data
+    GpuVector<double> result_vec(n);
+
+    // Test transfer AMGX -> GPU
+    AmgxInterface::updateGpuVectorFromAmgx(amgx_vec, result_vec);
+
+    // Verify data integrity
+    std::vector<double> result_data = result_vec.asStdVector();
+    BOOST_CHECK_EQUAL(result_data.size(), host_data.size());
+    for (size_t i = 0; i < n; ++i) {
+        BOOST_CHECK_CLOSE(result_data[i], host_data[i], 1e-10);
+    }
+
+    // Clean up
+    AmgxInterface::destroyVector(amgx_vec);
+}
+
+// Test matrix transfer functionality
+BOOST_AUTO_TEST_CASE(TestMatrixTransfer)
+{
+    // Setup AMGX resources
+    AmgxResourceGuard resources;
+
+    // Create a simple tridiagonal matrix
+    const int N = 10;
+    const int blockSize = 1;
+
+    using FieldMatrix = Dune::FieldMatrix<double, blockSize, blockSize>;
+    using BCRSMatrix = Dune::BCRSMatrix<FieldMatrix>;
+
+    BCRSMatrix dune_matrix(N, N, BCRSMatrix::row_wise);
+
+    // Set up matrix structure (tridiagonal)
+    for (auto row = dune_matrix.createbegin(); row != dune_matrix.createend(); ++row) {
+        const int i = row.index();
+        if (i > 0) {
+            row.insert(i - 1); // Lower diagonal
+        }
+        row.insert(i); // Diagonal
+        if (i < N - 1) {
+            row.insert(i + 1); // Upper diagonal
+        }
+    }
+
+    // Fill with values
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<double> dist(-10.0, 10.0);
+
+    for (int i = 0; i < N; ++i) {
+        dune_matrix[i][i] = 10.0 + std::abs(dist(gen)); // Diagonally dominant
+
+        if (i > 0) {
+            dune_matrix[i][i - 1] = dist(gen);
+        }
+
+        if (i < N - 1) {
+            dune_matrix[i][i + 1] = dist(gen);
+        }
+    }
+
+    // Convert to GpuSparseMatrix
+    auto gpu_matrix = GpuSparseMatrix<double>::fromMatrix(dune_matrix);
+
+    // Create AMGX matrix and upload data
+    AMGX_matrix_handle amgx_matrix = AmgxInterface::createMatrix(resources.getResources(), AMGX_mode_dDDI);
+    AmgxInterface::updateAmgxMatrixFromGpuSparseMatrix(gpu_matrix, amgx_matrix);
+
+    // Verify matrix dimensions after upload
+    int n, block_dimx, block_dimy;
+    OPM_AMGX_SAFE_CALL(AMGX_matrix_get_size(amgx_matrix, &n, &block_dimx, &block_dimy));
+    BOOST_CHECK_EQUAL(n, N);
+    BOOST_CHECK_EQUAL(block_dimx, blockSize);
+    BOOST_CHECK_EQUAL(block_dimy, blockSize);
+
+    // Test round-trip: upload to AMGX, then download back
+    auto gpu_matrix_copy = GpuSparseMatrix<double>::fromMatrix(dune_matrix);
+    AmgxInterface::updateGpuSparseMatrixFromAmgxMatrix(amgx_matrix, gpu_matrix_copy);
+
+    // Verify round-trip preserves original values
+    const auto original_values = gpu_matrix.getNonZeroValues().asStdVector();
+    const auto downloaded_values = gpu_matrix_copy.getNonZeroValues().asStdVector();
+    BOOST_CHECK_EQUAL(downloaded_values.size(), original_values.size());
+    for (std::size_t i = 0; i < original_values.size(); ++i) {
+        BOOST_CHECK_CLOSE(downloaded_values[i], original_values[i], 1e-10);
+    }
+
+    // Test coefficient update: modify GPU matrix and update AMGX coefficients only
+    gpu_matrix.getNonZeroValues() *= 2.0;
+    AmgxInterface::updateAmgxMatrixCoefficientsFromGpuSparseMatrix(gpu_matrix, amgx_matrix);
+
+    // Verify coefficient update worked
+    AmgxInterface::updateGpuSparseMatrixFromAmgxMatrix(amgx_matrix, gpu_matrix_copy);
+    const auto updated_values = gpu_matrix_copy.getNonZeroValues().asStdVector();
+    for (std::size_t i = 0; i < original_values.size(); ++i) {
+        BOOST_CHECK_CLOSE(updated_values[i], 2.0 * original_values[i], 1e-10);
+    }
+
+    // Clean up
+    AmgxInterface::destroyMatrix(amgx_matrix, gpu_matrix);
+}
+
+// Test templated utility functions
+BOOST_AUTO_TEST_CASE(TestTemplatedUtilities)
+{
+    // Test is_gpu_type trait
+    BOOST_CHECK(is_gpu_type<GpuVector<double>>::value);
+    BOOST_CHECK(is_gpu_type<GpuSparseMatrix<double>>::value);
+    BOOST_CHECK(!is_gpu_type<std::vector<double>>::value);
+
+    // Test AMGX mode determination
+    AMGX_Mode mode_dd = AmgxInterface::determineAmgxMode<double, double>();
+    BOOST_CHECK_EQUAL(mode_dd, AMGX_mode_dDDI);
+
+    AMGX_Mode mode_ff = AmgxInterface::determineAmgxMode<float, float>();
+    BOOST_CHECK_EQUAL(mode_ff, AMGX_mode_dFFI);
+}
+
+// Test error handling
+BOOST_AUTO_TEST_CASE(TestErrorHandling)
+{
+    // Setup AMGX resources
+    AmgxResourceGuard resources;
+
+    // Create vectors of different sizes
+    GpuVector<double> vec1(10);
+    GpuVector<double> vec2(5);
+
+    // Create AMGX vector
+    AMGX_vector_handle amgx_vec
+        = AmgxInterface::createVector(resources.getResources(), AmgxInterface::determineAmgxMode<double, double>());
+
+    // Set dimensions for the AMGX vector
+    OPM_AMGX_SAFE_CALL(AMGX_vector_set_zero(amgx_vec, 10, 1));
+
+    // This should work fine
+    AmgxInterface::updateAmgxFromGpuVector(vec1, amgx_vec);
+
+    // This should throw an exception due to size mismatch
+    BOOST_CHECK_THROW(AmgxInterface::updateGpuVectorFromAmgx(amgx_vec, vec2), AmgxError);
+
+    // Clean up
+    AmgxInterface::destroyVector(amgx_vec);
+}
+
+// Custom main function to properly initialize/finalize AMGX
+int
+main(int argc, char** argv)
+{
+    return boost::unit_test::unit_test_main(init_unit_test, argc, argv);
+}

tests/test_AmgxPreconditioner.cpp

@@ -29,30 +29,56 @@
 #include "AmgxPreconditionerTestHelper.hpp"
 
 
-BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerMatDoubleVecDouble)
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerCpuInputFloat)
 {
-    testAmgxPreconditioner<double, double>();
+    testAmgxPreconditioner<float, false>();
 }
 
-// This test is disabled because it crashes the program with the following error:
-// "Mixed precision modes not currently supported for CUDA 10.1 or later."
-//BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerMatFloatVecDouble)
-//{
-//    testAmgxPreconditioner<float, double>();
-//}
-
-
-BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerMatFloatVecFloat)
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerCpuInputDouble)
 {
-    testAmgxPreconditioner<float, float>();
+    testAmgxPreconditioner<double, false>();
 }
 
-bool init_unit_test_func()
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerUpdateCpuInputFloat)
+{
+    testAmgxPreconditionerUpdate<float, false>();
+}
+
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerUpdateCpuInputDouble)
+{
+    testAmgxPreconditionerUpdate<double, false>();
+}
+
+
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerGpuInputFloat)
+{
+    testAmgxPreconditioner<float, true>();
+}
+
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerGpuInputDouble)
+{
+    testAmgxPreconditioner<double, true>();
+}
+
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerUpdateGpuInputFloat)
+{
+    testAmgxPreconditionerUpdate<float, true>();
+}
+
+BOOST_AUTO_TEST_CASE(TestAmgxPreconditionerUpdateGpuInputDouble)
+{
+    testAmgxPreconditionerUpdate<double, true>();
+}
+
+
+bool
+init_unit_test_func()
 {
     return true;
 }
 
-int main(int argc, char** argv)
+int
+main(int argc, char** argv)
 {
     AMGX_SAFE_CALL(AMGX_initialize());