diff --git a/.github/workflows/continuous.yml b/.github/workflows/continuous.yml index 5f5f524..71ceafd 100644 --- a/.github/workflows/continuous.yml +++ b/.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" diff --git a/CMakeLists.txt b/CMakeLists.txt index 474ef9b..aad2499 100755 --- a/CMakeLists.txt +++ b/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) diff --git a/cmake/recipes/trilinos.cmake b/cmake/recipes/trilinos.cmake new file mode 100644 index 0000000..15174ea --- /dev/null +++ b/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() \ No newline at end of file diff --git a/linear-solver-spec.json b/linear-solver-spec.json index 1e59528..330aa0f 100644 --- a/linear-solver-spec.json +++ b/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." } ] \ No newline at end of file diff --git a/src/polysolve/linear/CMakeLists.txt b/src/polysolve/linear/CMakeLists.txt index a742a35..47d60a5 100644 --- a/src/polysolve/linear/CMakeLists.txt +++ b/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}) diff --git a/src/polysolve/linear/Solver.cpp b/src/polysolve/linear/Solver.cpp index 25595be..2ea1921 100644 --- a/src/polysolve/linear/Solver.cpp +++ b/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(); #endif +#ifdef POLYSOLVE_WITH_TRILINOS + } + else if (solver == "Trilinos") + { + return std::make_unique(); +#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", diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index 58e281f..3df9128 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -4,6 +4,7 @@ #include + #define POLYSOLVE_DELETE_MOVE_COPY(Base) \ Base(Base &&) = delete; \ Base &operator=(Base &&) = delete; \ @@ -34,6 +35,7 @@ namespace polysolve::linear public: // Shortcut alias typedef Eigen::VectorXd VectorXd; + typedef Eigen::MatrixXd MatrixXd; template using Ref = Eigen::Ref; @@ -123,6 +125,7 @@ namespace polysolve::linear /// and initialized. } /// virtual void solve(const Ref b, Ref x) = 0; + virtual void solve(const Ref b, const Ref nullspace, Ref x) {} /// @brief Name of the solver type (for debugging purposes) virtual std::string name() const { return ""; } diff --git a/src/polysolve/linear/TrilinosSolver.cpp b/src/polysolve/linear/TrilinosSolver.cpp new file mode 100644 index 0000000..4e94306 --- /dev/null +++ b/src/polysolve/linear/TrilinosSolver.cpp @@ -0,0 +1,288 @@ +#ifdef POLYSOLVE_WITH_TRILINOS + +//////////////////////////////////////////////////////////////////////////////// +#include "TrilinosSolver.hpp" +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace polysolve::linear +{ + TrilinosSolver::TrilinosSolver() + { + precond_num_ = 0; + int done_already; + MPI_Initialized(&done_already); + if (!done_already) + { + /* Initialize MPI */ + int argc = 1; + char name[] = ""; + char *argv[] = {name}; + char **argvv = &argv[0]; + MPI_Init(&argc, &argvv); + CommPtr = new Epetra_MpiComm(MPI_COMM_WORLD); + } + } + + //////////////////////////////////////////////////////////////// + void TrilinosSolver::set_parameters(const json ¶ms) + { + if (params.contains("Trilinos")) + { + if (params["Trilinos"].contains("block_size")) + { + if (params["Trilinos"]["block_size"] == 2 || params["Trilinos"]["block_size"] == 3) + { + numPDEs = params["Trilinos"]["block_size"]; + } + } + if (params["Trilinos"].contains("max_iter")) + { + max_iter_ = params["Trilinos"]["max_iter"]; + } + if (params["Trilinos"].contains("tolerance")) + { + conv_tol_ = params["Trilinos"]["tolerance"]; + } + } + } + + ///////////////////////////////////////////////// + void TrilinosSolver::get_info(json ¶ms) const + { + params["num_iterations"] = iterations_; + params["final_res_norm"] = residual_error_; + } + + ///////////////////////////////////////////////// + void TrilinosSolver::factorize(const StiffnessMatrix &Ain) + { + assert(precond_num_ > 0); + // Eigen::saveMarket(Ain,"/home/yiwei/matrix_struct/A_nonLinear.mtx"); + // Eigen::saveMarket(test_vertices,"/home/yiwei/matrix_struct/vec.mtx"); + Eigen::SparseMatrix Arow(Ain); + int mypid = CommPtr->MyPID(); + int indexBase=0; + int numGlobalElements = Arow.nonZeros(); + int numGlobalRows=Arow.rows(); + + int numNodes= numGlobalRows /numPDEs; + if ((numGlobalRows - numNodes * numPDEs) != 0 && !mypid){ + throw std::runtime_error("Number of matrix rows is not divisible by #dofs"); + } + int numMyNodes; + int nproc = CommPtr->NumProc(); + if (CommPtr->MyPID() < nproc-1) numMyNodes = numNodes / nproc; + else numMyNodes = numNodes - (numNodes/nproc) * (nproc-1); + delete rowMap; + delete A; + rowMap = new Epetra_Map(numGlobalRows,numMyNodes*numPDEs,indexBase,(*CommPtr)); + + A = new Epetra_CrsMatrix(Copy,*rowMap,0); //Can allocate memory for each row in advance + + { + int nnzs=0; + for (int k=0 ; k < Arow.outerSize(); k++) + { + // std::cout<InsertGlobalValues(k,numEntries,values,indices); + nnzs=nnzs+numEntries; + } + } + A->FillComplete(); + } + + namespace + { + void TrilinosML_SetDefaultOptions(Teuchos::ParameterList &MLList) + { + std::string aggr_type="Uncoupled-MIS"; + double str_connect=0.08; + ML_Epetra::SetDefaults("SA", MLList); + MLList.set("aggregation: type",aggr_type); // Aggresive Method + // MLList.set("aggregation: type","Uncoupled"); // Fixed size + + + MLList.set("aggregation: threshold",str_connect); + + // Smoother Settings + MLList.set("smoother: type","Chebyshev"); + MLList.set("smoother: sweeps", 5); //Chebyshev degree + MLList.set("smoother: Chebyshev alpha",30.0); + + //Coarser Settings + MLList.set("coarse: max size",1000); + + MLList.set("ML output",0); + } + } + + + void TrilinosSolver::solve(const Eigen::Ref rhs, const Eigen::Ref nullspace, Eigen::Ref result) + { + int output=10; //how often to print residual history + Teuchos::ParameterList MLList; + TrilinosML_SetDefaultOptions(MLList); + MLList.set("PDE equations",numPDEs); + + if (nullspace.cols()==3) + { + Eigen::VectorXd col1 = nullspace.col(0); + int size1 = col1.size(); + double* arr1 = new double[size1]; + for(int i = 0; i < size1; ++i) { + arr1[i] = col1(i); + } + Eigen::VectorXd col2 = nullspace.col(1); + int size2 = col2.size(); + double* arr2 = new double[size2]; + for(int i = 0; i < size2; ++i) { + arr2[i] = col2(i); + } + Eigen::VectorXd col3 = nullspace.col(2); + int size3 = col3.size(); + double* arr3 = new double[size3]; + for(int i = 0; i < size3; ++i) { + arr3[i] = col3(i); + } + MLList.set("null space: type","elasticity from coordinates"); + MLList.set("x-coordinates", arr1); // nullspace.col(0).data() + MLList.set("y-coordinates", arr2); // nullspace.col(1).data() + MLList.set("z-coordinates", arr3); // nullspace.col(2).data() + MLList.set("aggregation: threshold",0.00); + } + if (nullspace.cols()==2) + { + Eigen::VectorXd col1 = nullspace.col(0); + int size1 = col1.size(); + double* arr1 = new double[size1]; + for(int i = 0; i < size1; ++i) { + arr1[i] = col1(i); + } + Eigen::VectorXd col2 = nullspace.col(1); + int size2 = col2.size(); + double* arr2 = new double[size2]; + for(int i = 0; i < size2; ++i) { + arr2[i] = col2(i); + } + MLList.set("null space: type","elasticity from coordinates"); + MLList.set("x-coordinates", arr1); // nullspace.col(0).data() + MLList.set("y-coordinates", arr2); // nullspace.col(1).data() + MLList.set("aggregation: threshold",0.00); + } + + delete MLPrec; + MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList); + Epetra_Vector x(A->RowMap()); + Epetra_Vector b(A->RowMap()); + for (size_t i = 0; i < rhs.size(); i++) + { + x[i]=result[i]; + b[i]=rhs[i]; + } + + Epetra_LinearProblem Problem(A,&x,&b); + AztecOO solver(Problem); + solver.SetAztecOption(AZ_solver, AZ_cg); + solver.SetPrecOperator(MLPrec); + solver.SetAztecOption(AZ_output, AZ_last); + + int status= solver.Iterate(max_iter_, conv_tol_ ); + if (status!=0 && status!=-3) + { + throw std::runtime_error("Early termination, not SPD"); + } + + if (CommPtr->MyPID() == 0) + { + residual_error_=solver.ScaledResidual (); + iterations_=solver.NumIters(); + } + + for (size_t i = 0; i < rhs.size(); i++) + { + result[i]=x[i]; + } + + } + + void TrilinosSolver::solve(const Eigen::Ref rhs, Eigen::Ref result) + { + int output=10; //how often to print residual history + Teuchos::ParameterList MLList; + TrilinosML_SetDefaultOptions(MLList); + MLList.set("PDE equations",numPDEs); + + delete MLPrec; + MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList); + Epetra_Vector x(A->RowMap()); + Epetra_Vector b(A->RowMap()); + for (size_t i = 0; i < rhs.size(); i++) + { + x[i]=result[i]; + b[i]=rhs[i]; + } + // std::cout<<"x[0] "<RowMap()); + // double mynorm; + // A->Multiply(false,x,workvec); + // workvec.Update(1.0,b,-1.0); + // b.Norm2(&mynorm); + // workvec.Scale(1./mynorm); + // workvec.Norm2(&mynorm); + if (CommPtr->MyPID() == 0) + { + + // std::cout<<"Max iterations "<175) + // { + // exit(); + // } + + + + for (size_t i = 0; i < rhs.size(); i++) + { + result[i]=x[i]; + } + + } + + TrilinosSolver:: ~TrilinosSolver() + { + delete A; + delete rowMap; + delete MLPrec; + MPI_Finalize() ; + } +} + +#endif diff --git a/src/polysolve/linear/TrilinosSolver.hpp b/src/polysolve/linear/TrilinosSolver.hpp new file mode 100644 index 0000000..0b3fdb2 --- /dev/null +++ b/src/polysolve/linear/TrilinosSolver.hpp @@ -0,0 +1,102 @@ +#pragma once + +#ifdef POLYSOLVE_WITH_TRILINOS + +//////////////////////////////////////////////////////////////////////////////// +#include "Solver.hpp" + +#include +#include +#include + +#ifdef HAVE_MPI +#include "mpi.h" +#include "Epetra_MpiComm.h" +#else +#include "Epetra_SerialComm.h" +#endif +#include "Teuchos_CommandLineProcessor.hpp" +#include "Epetra_Map.h" +#include "Epetra_Vector.h" +#include "Epetra_CrsMatrix.h" +#include "Epetra_LinearProblem.h" +#include "EpetraExt_BlockMapIn.h" +#include "EpetraExt_CrsMatrixIn.h" +#include "EpetraExt_RowMatrixOut.h" +#include "EpetraExt_MultiVectorOut.h" +#include "EpetraExt_MultiVectorIn.h" +#include "AztecOO.h" + +#include "ml_include.h" +#include "ml_MultiLevelPreconditioner.h" +#include "ml_epetra.h" +#include + +//////////////////////////////////////////////////////////////////////////////// +// +// WARNING: +// The matrix is assumed to be in row-major format, since AMGCL assumes that the +// outer index is for row. If the matrix is symmetric, you are fine, because CSR +// and CSC are the same. If the matrix is not symmetric and you pass in a +// column-major matrix, the solver will actually solve A^T x = b. +// + +namespace polysolve::linear +{ + + class TrilinosSolver : public Solver + { + + public: + TrilinosSolver(); + ~TrilinosSolver(); + + private: + POLYSOLVE_DELETE_MOVE_COPY(TrilinosSolver) + + public: + ////////////////////// + // Public interface // + ////////////////////// + + // Set solver parameters + virtual void set_parameters(const json ¶ms) override; + + // Retrieve memory information from Pardiso + virtual void get_info(json ¶ms) const override; + + // Analyze sparsity pattern + virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) override { precond_num_ = precond_num; } + + // Factorize system matrix + virtual void factorize(const StiffnessMatrix &A) override; + + // Solve the linear system Ax = b + virtual void solve(const Ref b, Ref x) override; + virtual void solve(const Ref b, const Ref nullspace, Ref x); + + // Name of the solver type (for debugging purposes) + virtual std::string name() const override { return "Trilinos AztecOO and ML"; } + + protected: + int numPDEs = 1; // 1 = scalar (Laplace), 2 or 3 = vector (Elasticity) + int max_iter_ = 1000; + double conv_tol_ = 1e-8; + size_t iterations_; + double residual_error_; + ML_Epetra::MultiLevelPreconditioner* MLPrec=NULL; + Epetra_Map *rowMap=NULL; + + private: + int precond_num_; + Epetra_CrsMatrix *A=NULL; +#ifdef HAVE_MPI + Epetra_MpiComm *CommPtr; +#else + Epetra_SerialComm *CommPtr; +#endif + }; + +} // namespace polysolve + +#endif \ No newline at end of file