Add Trilinos for linear solver

.github/workflows/continuous.yml

@@ -36,12 +36,13 @@ jobs:
       - name: Dependencies (Linux)
         if: runner.os == 'Linux'
         run: |
-          sudo apt update
-          sudo apt -o Acquire::Retries=3 install \
+          sudo apt-get update
+          sudo apt-get -o Acquire::Retries=3 install \
           libblas-dev \
           libglu1-mesa-dev \
           xorg-dev \
           mpi \
+          trilinos-dev \
           ccache
           echo 'CACHE_PATH=~/.ccache' >> "$GITHUB_ENV"
 

CMakeLists.txt

@@ -77,9 +77,11 @@ option(POLYSOLVE_WITH_SUPERLU       "Enable SuperLU library"
 option(POLYSOLVE_WITH_MKL           "Enable MKL library"  ${POLYSOLVE_NOT_ON_APPLE_SILICON})
 option(POLYSOLVE_WITH_CUSOLVER      "Enable cuSOLVER library"                           OFF)
 option(POLYSOLVE_WITH_PARDISO       "Enable Pardiso library"                            OFF)
-option(POLYSOLVE_WITH_HYPRE         "Enable hypre"                                       ON)
+option(POLYSOLVE_WITH_HYPRE         "Enable hypre"                                      OFF)
 option(POLYSOLVE_WITH_AMGCL         "Use AMGCL"                                          ON)
 option(POLYSOLVE_WITH_SPECTRA       "Enable Spectra library"                             ON)
+option(POLYSOLVE_WITH_TRILINOS      "Enable Trilinos"                                    ON)
+option(BUILD_TRILINOS_FROM_SOURCE   "Build trilinos from source"                        OFF)
 
 # Sanitizer options
 option(POLYSOLVE_SANITIZE_ADDRESS   "Sanitize Address"                                  OFF)
@@ -217,7 +219,7 @@ include(jse)
 target_link_libraries(polysolve_linear PUBLIC jse::jse)
 
 # Hypre (GNU Lesser General Public License)
-if(POLYSOLVE_WITH_HYPRE)
+if(POLYSOLVE_WITH_HYPRE AND NOT POLYSOLVE_WITH_TRILINOS)
     include(hypre)
     target_link_libraries(polysolve_linear PUBLIC HYPRE::HYPRE)
     target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_WITH_HYPRE)
@@ -283,6 +285,17 @@ if(POLYSOLVE_WITH_SPECTRA)
     target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_WITH_SPECTRA)
 endif()
 
+# Trilinos
+if (POLYSOLVE_WITH_TRILINOS)
+    include(trilinos)
+    if(TARGET Trilinos::Trilinos)
+        target_link_libraries(polysolve_linear PUBLIC Trilinos::Trilinos)
+        target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_TRILINOS)
+    else()
+        message(WARNING "Trilinos not found, solver will not be available.")
+    endif()
+endif()
+
 # cuSolver solvers
 if(POLYSOLVE_WITH_CUSOLVER)
     include(cusolverdn)

cmake/recipes/trilinos.cmake

@@ -0,0 +1,43 @@
+if(TARGET Trilinos::Trilinos)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'Trilinos::Trilinos'")
+
+find_package(Trilinos COMPONENTS ML Epetra)
+
+if(NOT Trilinos_FOUND)
+    set(POLYSOLVE_WITH_TRILINOS OFF)
+    message("Trilinos not found.")
+endif()
+
+find_package(MPI)
+
+if(NOT MPI_FOUND)
+    set(POLYSOLVE_WITH_TRILINOS OFF)
+    message("MPI not found.")
+endif()
+
+if(Trilinos_FOUND)
+    if(MPI_FOUND)
+        MESSAGE("\nFound Trilinos!  Here are the details: ")
+        MESSAGE("   Trilinos_DIR = ${Trilinos_DIR}")
+        MESSAGE("   Trilinos_VERSION = ${Trilinos_VERSION}")
+        MESSAGE("   Trilinos_PACKAGE_LIST = ${Trilinos_PACKAGE_LIST}")
+        MESSAGE("   Trilinos_LIBRARIES = ${Trilinos_LIBRARIES} ")
+        MESSAGE("   Trilinos_INCLUDE_DIRS = ${Trilinos_INCLUDE_DIRS} ")
+        MESSAGE("   Trilinos_TPL_LIST = ${Trilinos_TPL_LIST}")
+        MESSAGE("   Trilinos_TPL_LIBRARIES = ${Trilinos_TPL_LIBRARIES}")
+        MESSAGE("   Trilinos_BUILD_SHARED_LIBS = ${Trilinos_BUILD_SHARED_LIBS}")
+        MESSAGE("End of Trilinos details\n")
+        # include(trilinos)
+        if(TARGET Trilinos::Trilinos)
+        else()
+            add_library(trilinos INTERFACE)
+            add_library(Trilinos::Trilinos ALIAS trilinos)
+            target_include_directories(trilinos INTERFACE ${Trilinos_INCLUDE_DIRS} )
+            target_link_libraries(trilinos INTERFACE ${Trilinos_LIBRARIES} )
+            target_link_libraries(trilinos INTERFACE MPI::MPI_C MPI::MPI_CXX )
+        endif()
+    endif()
+endif()

linear-solver-spec.json

@@ -15,7 +15,8 @@
             "Eigen::MINRES",
             "Pardiso",
             "Hypre",
-            "AMGCL"
+            "AMGCL",
+            "Trilinos"
         ],
         "doc": "Settings for the linear solver."
     },
@@ -42,6 +43,7 @@
             "Pardiso",
             "Hypre",
             "AMGCL",
+            "Trilinos",
             "Eigen::LeastSquaresConjugateGradient",
             "Eigen::DGMRES",
             "Eigen::ConjugateGradient",
@@ -153,6 +155,17 @@
         ],
         "doc": "Settings for the AMGCL solver."
     },
+    {
+        "pointer": "/Trilinos",
+        "default": null,
+        "type": "object",
+        "optional": [
+            "max_iter",
+            "tolerance",
+            "block_size"
+        ],
+        "doc": "Settings for the Trilinos solver."
+    },
     {
         "pointer": "/Eigen::LeastSquaresConjugateGradient/max_iter",
         "default": 1000,
@@ -421,5 +434,29 @@
         "default": 0,
         "type": "float",
         "doc": "Aggregation epsilon strong."
+    },
+    {
+        "pointer": "/Trilinos/max_iter",
+        "default": 1000,
+        "type": "int",
+        "doc": "Maximum number of iterations."
+    },
+    {
+        "pointer": "/Trilinos/block_size",
+        "default": 3,
+        "type": "int",
+        "doc": "Aggregation epsilon strong."
+    },
+    {
+        "pointer": "/Trilinos/tolerance",
+        "default": 1e-8,
+        "type": "float",
+        "doc": "Convergence tolerance."
+    },
+    {
+        "pointer": "/Trilinos/is_nullspace",
+        "default": false,
+        "type": "bool",
+        "doc": "Is nullspace or not."
     }
 ]

src/polysolve/linear/CMakeLists.txt

@@ -15,6 +15,8 @@ set(SOURCES
     Pardiso.hpp
     SaddlePointSolver.cpp
     SaddlePointSolver.hpp
+    TrilinosSolver.cpp
+    TrilinosSolver.hpp
 )
 
 source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES})

src/polysolve/linear/Solver.cpp

@@ -36,6 +36,9 @@
 #ifdef POLYSOLVE_WITH_AMGCL
 #include "AMGCL.hpp"
 #endif
+#ifdef POLYSOLVE_WITH_TRILINOS
+#include "TrilinosSolver.hpp"
+#endif
 #ifdef POLYSOLVE_WITH_CUSOLVER
 #include "CuSolverDN.cuh"
 #endif
@@ -377,6 +380,12 @@ namespace polysolve::linear
         {
             return std::make_unique<AMGCL>();
 #endif
+#ifdef POLYSOLVE_WITH_TRILINOS
+        }
+        else if (solver == "Trilinos")
+        {
+            return std::make_unique<TrilinosSolver>();
+#endif
 #if EIGEN_VERSION_AT_LEAST(3, 3, 0)
             // Available only with Eigen 3.3.0 and newer
 #ifndef POLYSOLVE_LARGE_INDEX
@@ -499,6 +508,9 @@ namespace polysolve::linear
 #ifdef POLYSOLVE_WITH_AMGCL
             "AMGCL",
 #endif
+#ifdef POLYSOLVE_WITH_TRILINOS
+            "Trilinos",
+#endif
 #if EIGEN_VERSION_AT_LEAST(3, 3, 0)
 #ifndef POLYSOLVE_LARGE_INDEX
             "Eigen::LeastSquaresConjugateGradient",

src/polysolve/linear/Solver.hpp

@@ -4,6 +4,7 @@
 
 #include <memory>
 
+
 #define POLYSOLVE_DELETE_MOVE_COPY(Base) \
     Base(Base &&) = delete;              \
     Base &operator=(Base &&) = delete;   \
@@ -34,6 +35,7 @@
     public:
         // Shortcut alias
         typedef Eigen::VectorXd VectorXd;
+        typedef Eigen::MatrixXd MatrixXd;
         template <typename T>
         using Ref = Eigen::Ref<T>;
 
@@ -123,6 +125,7 @@
         ///                      and initialized. }
         ///
         virtual void solve(const Ref<const VectorXd> b, Ref<VectorXd> x) = 0;
+        virtual void solve(const Ref<const VectorXd> b, const Ref<const MatrixXd> nullspace, Ref<VectorXd> x) {}
 
         /// @brief Name of the solver type (for debugging purposes)
         virtual std::string name() const { return ""; }