Add Trilinos for linear solver

CMakeLists.txt

@@ -80,6 +80,8 @@ option(POLYSOLVE_WITH_PARDISO       "Enable Pardiso library"
 option(POLYSOLVE_WITH_HYPRE         "Enable hypre"                                       ON)
 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,13 @@ if(POLYSOLVE_WITH_SPECTRA)
     target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_WITH_SPECTRA)
 endif()
 
+# Trilinos
+if (POLYSOLVE_WITH_TRILINOS)
+    include(trilinos)
+    target_link_libraries(polysolve_linear PUBLIC Trilinos::Trilinos)
+    target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_WITH_TRILINOS)
+endif()
+
 # cuSolver solvers
 if(POLYSOLVE_WITH_CUSOLVER)
     include(cusolverdn)

cmake/recipes/trilinos.cmake

@@ -0,0 +1,37 @@
+if(TARGET Trilinos::Trilinos)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'Trilinos::Trilinos'")
+
+find_package(Trilinos COMPONENTS ML Epetra)
+
+if(NOT Trilinos_FOUND)
+    message("Trilinos not found.")
+endif()
+
+find_package(MPI)
+
+if(NOT MPI_FOUND)
+    message("MPI not found.")
+endif()
+
+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()

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,16 @@
         ],
         "doc": "Settings for the AMGCL solver."
     },
+    {
+        "pointer": "/Trilinos",
+        "default": null,
+        "type": "object",
+        "optional": [
+            "solver",
+            "precond"
+        ],
+        "doc": "Settings for the Trilinos solver."
+    },
     {
         "pointer": "/Eigen::LeastSquaresConjugateGradient/max_iter",
         "default": 1000,

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,12 +36,42 @@
 #ifdef POLYSOLVE_WITH_AMGCL
 #include "AMGCL.hpp"
 #endif
+#ifdef POLYSOLVE_WITH_TRILINOS
+#include "TrilinosSolver.hpp"
+#endif
 #ifdef POLYSOLVE_WITH_CUSOLVER
 #include "CuSolverDN.cuh"
 #endif
 #include <unsupported/Eigen/IterativeSolvers>
 
 ////////////////////////////////////////////////////////////////////////////////
+Eigen::MatrixXd test_vertices;
+Eigen::MatrixXd init_vertices;
+std::vector<int> test_boundary_nodes;
+Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector<int> &boundary_nodes)
+{
+    // Remove boundary vertices
+    if (boundary_nodes.empty())
+    {
+        return vertices;
+    }
+    else
+    {
+        std::vector<int> order_nodes = boundary_nodes;
+        std::sort(order_nodes.begin(), order_nodes.end());
+        Eigen::MatrixXd out_vertices;
+        std::vector<int> keep;
+        for (int i = 0; i < vertices.rows(); i++)
+        {
+            if (!std::binary_search(order_nodes.begin(), order_nodes.end(),i))
+            {
+                keep.push_back(i);
+            }
+        }
+        out_vertices = vertices(keep, Eigen::all);
+        return out_vertices;
+    }
+}
 
 namespace polysolve::linear
 {
@@ -377,6 +407,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 +535,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,12 +4,21 @@
 
 #include <memory>
 
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <vector>
+
 #define POLYSOLVE_DELETE_MOVE_COPY(Base) \
     Base(Base &&) = delete;              \
     Base &operator=(Base &&) = delete;   \
     Base(const Base &) = delete;         \
     Base &operator=(const Base &) = delete;
 
+extern Eigen::MatrixXd init_vertices;
+extern Eigen::MatrixXd test_vertices;
+extern std::vector<int> test_boundary_nodes;
+Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector<int> &boundary_nodes);
+
 ////////////////////////////////////////////////////////////////////////////////
 // TODO:
 // - [ ] Support both RowMajor + ColumnMajor sparse matrices
@@ -25,6 +34,11 @@ namespace spdlog
 
 namespace polysolve::linear
 {
+#ifdef POLYSOLVE_LARGE_INDEX
+    typedef Eigen::SparseMatrix<double, Eigen::ColMajor, std::ptrdiff_t> StiffnessMatrix;
+#else
+    typedef Eigen::SparseMatrix<double, Eigen::ColMajor> StiffnessMatrix;
+#endif
     /**
      * @brief      Base class for linear solver.
      */