From a376e0d5726e802107cdb0f801d818ef5f5ce746 Mon Sep 17 00:00:00 2001 From: xDarkLemon Date: Wed, 14 Feb 2024 17:51:27 -0500 Subject: [PATCH 1/5] Add Trilinos --- CMakeLists.txt | 11 +- cmake/recipes/trilinos.cmake | 37 +++ linear-solver-spec.json | 14 +- src/polysolve/linear/CMakeLists.txt | 2 + src/polysolve/linear/Solver.cpp | 39 +++ src/polysolve/linear/Solver.hpp | 14 + src/polysolve/linear/TrilinosSolver.cpp | 351 ++++++++++++++++++++++++ src/polysolve/linear/TrilinosSolver.hpp | 103 +++++++ 8 files changed, 569 insertions(+), 2 deletions(-) create mode 100644 cmake/recipes/trilinos.cmake create mode 100644 src/polysolve/linear/TrilinosSolver.cpp create mode 100644 src/polysolve/linear/TrilinosSolver.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 587f37e1..a7d7f1f3 100755 --- a/CMakeLists.txt +++ b/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 PRIVATE -DPOLYSOLVE_WITH_HYPRE) @@ -283,6 +285,13 @@ if(POLYSOLVE_WITH_SPECTRA) target_compile_definitions(polysolve_linear PRIVATE -DPOLYSOLVE_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) diff --git a/cmake/recipes/trilinos.cmake b/cmake/recipes/trilinos.cmake new file mode 100644 index 00000000..1b15890e --- /dev/null +++ b/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() \ No newline at end of file diff --git a/linear-solver-spec.json b/linear-solver-spec.json index 1e59528b..c335bd4a 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,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, diff --git a/src/polysolve/linear/CMakeLists.txt b/src/polysolve/linear/CMakeLists.txt index a742a357..47d60a5c 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 c1e8a4e4..761237a3 100644 --- a/src/polysolve/linear/Solver.cpp +++ b/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 //////////////////////////////////////////////////////////////////////////////// +Eigen::MatrixXd test_vertices; +Eigen::MatrixXd init_vertices; +std::vector test_boundary_nodes; +Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes) +{ + // Remove boundary vertices + if (boundary_nodes.empty()) + { + return vertices; + } + else + { + std::vector order_nodes = boundary_nodes; + std::sort(order_nodes.begin(), order_nodes.end()); + Eigen::MatrixXd out_vertices; + std::vector 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(); #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 +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", diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index deaae018..841250d7 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -4,12 +4,21 @@ #include +#include +#include +#include + #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 test_boundary_nodes; +Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &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 StiffnessMatrix; +#else + typedef Eigen::SparseMatrix StiffnessMatrix; +#endif /** * @brief Base class for linear solver. */ diff --git a/src/polysolve/linear/TrilinosSolver.cpp b/src/polysolve/linear/TrilinosSolver.cpp new file mode 100644 index 00000000..0f4e36c9 --- /dev/null +++ b/src/polysolve/linear/TrilinosSolver.cpp @@ -0,0 +1,351 @@ +#ifdef POLYSOLVE_WITH_TRILINOS + +//////////////////////////////////////////////////////////////////////////////// +#include "TrilinosSolver.hpp" +#include +#include +#include +/////////////////////////////////s/////////////////////////////////////////////// + +namespace polysolve::linear +{ + namespace{ + //////////////////////////////////////////////////////////////// + // int rigid_body_mode(int ndim, const std::vector &coo, std::vector &B, bool transpose = true) { + + // size_t n = coo.size(); + // int nmodes = (ndim == 2 ? 3 : 6); + // B.resize(n * nmodes, 0.0); + + // const int stride1 = transpose ? 1 : nmodes; + // const int stride2 = transpose ? n : 1; + // // int stride1 = nmodes; + // // int stride2 = 1; + + // double sn = 1 / sqrt(n); + + // if (ndim == 2) { + // for(size_t i = 0; i < n; ++i) { + // size_t nod = i / ndim; + // size_t dim = i % ndim; + + // double x = coo[nod * 2 + 0]; + // double y = coo[nod * 2 + 1]; + + // // Translation + // B[i * stride1 + dim * stride2] = sn; + + // // Rotation + // switch(dim) { + // case 0: + // B[i * stride1 + 2 * stride2] = -y; + // break; + // case 1: + // B[i * stride1 + 2 * stride2] = x; + // break; + // } + // } + // } else if (ndim == 3) { + // for(size_t i = 0; i < n; ++i) { + // size_t nod = i / ndim; + // size_t dim = i % ndim; + + // double x = coo[nod * 3 + 0]; + // double y = coo[nod * 3 + 1]; + // double z = coo[nod * 3 + 2]; + + // // Translation + // B[i * stride1 + dim * stride2] = sn; + + // // Rotation + // switch(dim) { + // case 0: + // B[i * stride1 + 5 * stride2] = -y; + // B[i * stride1 + 4 * stride2] = z; + // break; + // case 1: + // B[i * stride1 + 5 * stride2] = x; + // B[i * stride1 + 3 * stride2] = -z; + // break; + // case 2: + // B[i * stride1 + 3 * stride2] = y; + // B[i * stride1 + 4 * stride2] = -x; + // break; + // } + // } + // } + + // // Orthonormalization + // std::array dot; + // for(int i = ndim; i < nmodes; ++i) { + // std::fill(dot.begin(), dot.end(), 0.0); + // for(size_t j = 0; j < n; ++j) { + // for(int k = 0; k < i; ++k) + // dot[k] += B[j * stride1 + k * stride2] * B[j * stride1 + i * stride2]; + // } + // double s = 0.0; + // for(size_t j = 0; j < n; ++j) { + // for(int k = 0; k < i; ++k) + // B[j * stride1 + i * stride2] -= dot[k] * B[j * stride1 + k * stride2]; + // s += B[j * stride1 + i * stride2] * B[j * stride1 + i * stride2]; + // } + // s = sqrt(s); + // for(size_t j = 0; j < n; ++j) + // B[j * stride1 + i * stride2] /= s; + // } + // return nmodes; + // } + // Produce a vector of rigid body modes + + } + TrilinosSolver::TrilinosSolver() + { + precond_num_ = 0; +#ifdef HAVE_MPI + 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); + } +#else + CommPtr=new Epetra_SerialComm; +#endif + } + + //////////////////////////////////////////////////////////////// + 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"]; + } + if (params["Trilinos"].contains("is_nullspace")) + { + is_nullspace_ = params["Trilinos"]["is_nullspace"]; + } + } + } + + ///////////////////////////////////////////////// + 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, Eigen::Ref result) + { + int output=10; //how often to print residual history + Teuchos::ParameterList MLList; + TrilinosML_SetDefaultOptions(MLList); + MLList.set("PDE equations",numPDEs); + + //Set null space + + // if (true) + // { + // int n=test_vertices.rows(); + // int NRbm=0; + // int NscalarDof=0; + + // if (numPDEs==2) + // { + // NRbm=3; + // rbm=new double[n*(NRbm+NscalarDof)*(numPDEs+NscalarDof)]; + // std::vector z_coord(n,0); + // ML_Coord2RBM(n,test_vertices.col(0).data(),test_vertices.col(1).data(),z_coord.data(),rbm,numPDEs,NscalarDof); + // } + // else + // { + // NRbm=6; + // rbm=new double[n*(NRbm+NscalarDof)*(numPDEs+NscalarDof)]; + // ML_Coord2RBM(n,test_vertices.col(0).data(),test_vertices.col(1).data(),test_vertices.col(2).data(),rbm,numPDEs,NscalarDof); + // } + + // MLList.set("null space: vectors",rbm); + // MLList.set("null space: dimension", NRbm); + // MLList.set("null space: type", "pre-computed"); + // MLList.set("aggregation: threshold",0.00); + // } + + /////////////////////////////////////////////////////////////////////// + + if (is_nullspace_) + { + if (test_vertices.cols()==3) + { + reduced_vertices=remove_boundary_vertices(test_vertices,test_boundary_nodes); + MLList.set("null space: type","elasticity from coordinates"); + MLList.set("x-coordinates", reduced_vertices.col(0).data()); + MLList.set("y-coordinates", reduced_vertices.col(1).data()); + MLList.set("z-coordinates", reduced_vertices.col(2).data()); + MLList.set("aggregation: threshold",0.00); + } + if (test_vertices.cols()==2) + { + reduced_vertices=remove_boundary_vertices(test_vertices,test_boundary_nodes); + MLList.set("null space: type","elasticity from coordinates"); + MLList.set("x-coordinates", reduced_vertices.col(0).data()); + MLList.set("y-coordinates", reduced_vertices.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]; + } + // 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; +#ifdef HAVE_MPI + MPI_Finalize() ; +#endif + } +} + +#endif diff --git a/src/polysolve/linear/TrilinosSolver.hpp b/src/polysolve/linear/TrilinosSolver.hpp new file mode 100644 index 00000000..5385b668 --- /dev/null +++ b/src/polysolve/linear/TrilinosSolver.hpp @@ -0,0 +1,103 @@ +#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; + + // 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_; + bool is_nullspace_ = true; + Eigen::MatrixXd reduced_vertices; + 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 From 31539bb5753d8695b64bcfc49a4f6a156652621c Mon Sep 17 00:00:00 2001 From: xDarkLemon Date: Fri, 16 Feb 2024 14:13:05 -0500 Subject: [PATCH 2/5] Add Trilinos spec --- linear-solver-spec.json | 29 ++++- src/polysolve/linear/Solver.cpp | 27 ---- src/polysolve/linear/Solver.hpp | 10 -- src/polysolve/linear/TrilinosSolver.cpp | 161 ++++++------------------ src/polysolve/linear/TrilinosSolver.hpp | 5 + 5 files changed, 67 insertions(+), 165 deletions(-) diff --git a/linear-solver-spec.json b/linear-solver-spec.json index c335bd4a..330aa0fb 100644 --- a/linear-solver-spec.json +++ b/linear-solver-spec.json @@ -160,8 +160,9 @@ "default": null, "type": "object", "optional": [ - "solver", - "precond" + "max_iter", + "tolerance", + "block_size" ], "doc": "Settings for the Trilinos solver." }, @@ -433,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/Solver.cpp b/src/polysolve/linear/Solver.cpp index 2a320d39..65fe79bf 100644 --- a/src/polysolve/linear/Solver.cpp +++ b/src/polysolve/linear/Solver.cpp @@ -45,33 +45,6 @@ #include //////////////////////////////////////////////////////////////////////////////// -Eigen::MatrixXd test_vertices; -Eigen::MatrixXd init_vertices; -std::vector test_boundary_nodes; -Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes) -{ - // Remove boundary vertices - if (boundary_nodes.empty()) - { - return vertices; - } - else - { - std::vector order_nodes = boundary_nodes; - std::sort(order_nodes.begin(), order_nodes.end()); - Eigen::MatrixXd out_vertices; - std::vector 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 { diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index 841250d7..b2c6cf33 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -14,11 +14,6 @@ Base(const Base &) = delete; \ Base &operator=(const Base &) = delete; -extern Eigen::MatrixXd init_vertices; -extern Eigen::MatrixXd test_vertices; -extern std::vector test_boundary_nodes; -Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes); - //////////////////////////////////////////////////////////////////////////////// // TODO: // - [ ] Support both RowMajor + ColumnMajor sparse matrices @@ -34,11 +29,6 @@ namespace spdlog namespace polysolve::linear { -#ifdef POLYSOLVE_LARGE_INDEX - typedef Eigen::SparseMatrix StiffnessMatrix; -#else - typedef Eigen::SparseMatrix StiffnessMatrix; -#endif /** * @brief Base class for linear solver. */ diff --git a/src/polysolve/linear/TrilinosSolver.cpp b/src/polysolve/linear/TrilinosSolver.cpp index 0f4e36c9..4be30094 100644 --- a/src/polysolve/linear/TrilinosSolver.cpp +++ b/src/polysolve/linear/TrilinosSolver.cpp @@ -5,103 +5,42 @@ #include #include #include -/////////////////////////////////s/////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// -namespace polysolve::linear +Eigen::MatrixXd test_vertices; +Eigen::MatrixXd init_vertices; +std::vector test_boundary_nodes; +Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes) { - namespace{ - //////////////////////////////////////////////////////////////// - // int rigid_body_mode(int ndim, const std::vector &coo, std::vector &B, bool transpose = true) { - - // size_t n = coo.size(); - // int nmodes = (ndim == 2 ? 3 : 6); - // B.resize(n * nmodes, 0.0); - - // const int stride1 = transpose ? 1 : nmodes; - // const int stride2 = transpose ? n : 1; - // // int stride1 = nmodes; - // // int stride2 = 1; - - // double sn = 1 / sqrt(n); - - // if (ndim == 2) { - // for(size_t i = 0; i < n; ++i) { - // size_t nod = i / ndim; - // size_t dim = i % ndim; - - // double x = coo[nod * 2 + 0]; - // double y = coo[nod * 2 + 1]; - - // // Translation - // B[i * stride1 + dim * stride2] = sn; - - // // Rotation - // switch(dim) { - // case 0: - // B[i * stride1 + 2 * stride2] = -y; - // break; - // case 1: - // B[i * stride1 + 2 * stride2] = x; - // break; - // } - // } - // } else if (ndim == 3) { - // for(size_t i = 0; i < n; ++i) { - // size_t nod = i / ndim; - // size_t dim = i % ndim; - - // double x = coo[nod * 3 + 0]; - // double y = coo[nod * 3 + 1]; - // double z = coo[nod * 3 + 2]; - - // // Translation - // B[i * stride1 + dim * stride2] = sn; - - // // Rotation - // switch(dim) { - // case 0: - // B[i * stride1 + 5 * stride2] = -y; - // B[i * stride1 + 4 * stride2] = z; - // break; - // case 1: - // B[i * stride1 + 5 * stride2] = x; - // B[i * stride1 + 3 * stride2] = -z; - // break; - // case 2: - // B[i * stride1 + 3 * stride2] = y; - // B[i * stride1 + 4 * stride2] = -x; - // break; - // } - // } - // } - - // // Orthonormalization - // std::array dot; - // for(int i = ndim; i < nmodes; ++i) { - // std::fill(dot.begin(), dot.end(), 0.0); - // for(size_t j = 0; j < n; ++j) { - // for(int k = 0; k < i; ++k) - // dot[k] += B[j * stride1 + k * stride2] * B[j * stride1 + i * stride2]; - // } - // double s = 0.0; - // for(size_t j = 0; j < n; ++j) { - // for(int k = 0; k < i; ++k) - // B[j * stride1 + i * stride2] -= dot[k] * B[j * stride1 + k * stride2]; - // s += B[j * stride1 + i * stride2] * B[j * stride1 + i * stride2]; - // } - // s = sqrt(s); - // for(size_t j = 0; j < n; ++j) - // B[j * stride1 + i * stride2] /= s; - // } - // return nmodes; - // } - // Produce a vector of rigid body modes - + // Remove boundary vertices + if (boundary_nodes.empty()) + { + return vertices; } + else + { + std::vector order_nodes = boundary_nodes; + std::sort(order_nodes.begin(), order_nodes.end()); + Eigen::MatrixXd out_vertices; + std::vector 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 +{ TrilinosSolver::TrilinosSolver() { precond_num_ = 0; -#ifdef HAVE_MPI +// #ifdef HAVE_MPI int done_already; MPI_Initialized(&done_already); if (!done_already) @@ -114,9 +53,9 @@ namespace polysolve::linear MPI_Init(&argc, &argvv); CommPtr = new Epetra_MpiComm(MPI_COMM_WORLD); } -#else - CommPtr=new Epetra_SerialComm; -#endif +// #else + // CommPtr=new Epetra_SerialComm; +// #endif } //////////////////////////////////////////////////////////////// @@ -226,36 +165,6 @@ namespace polysolve::linear Teuchos::ParameterList MLList; TrilinosML_SetDefaultOptions(MLList); MLList.set("PDE equations",numPDEs); - - //Set null space - - // if (true) - // { - // int n=test_vertices.rows(); - // int NRbm=0; - // int NscalarDof=0; - - // if (numPDEs==2) - // { - // NRbm=3; - // rbm=new double[n*(NRbm+NscalarDof)*(numPDEs+NscalarDof)]; - // std::vector z_coord(n,0); - // ML_Coord2RBM(n,test_vertices.col(0).data(),test_vertices.col(1).data(),z_coord.data(),rbm,numPDEs,NscalarDof); - // } - // else - // { - // NRbm=6; - // rbm=new double[n*(NRbm+NscalarDof)*(numPDEs+NscalarDof)]; - // ML_Coord2RBM(n,test_vertices.col(0).data(),test_vertices.col(1).data(),test_vertices.col(2).data(),rbm,numPDEs,NscalarDof); - // } - - // MLList.set("null space: vectors",rbm); - // MLList.set("null space: dimension", NRbm); - // MLList.set("null space: type", "pre-computed"); - // MLList.set("aggregation: threshold",0.00); - // } - - /////////////////////////////////////////////////////////////////////// if (is_nullspace_) { @@ -342,9 +251,9 @@ namespace polysolve::linear delete A; delete rowMap; delete MLPrec; -#ifdef HAVE_MPI +// #ifdef HAVE_MPI MPI_Finalize() ; -#endif +// #endif } } diff --git a/src/polysolve/linear/TrilinosSolver.hpp b/src/polysolve/linear/TrilinosSolver.hpp index 5385b668..56d839c6 100644 --- a/src/polysolve/linear/TrilinosSolver.hpp +++ b/src/polysolve/linear/TrilinosSolver.hpp @@ -41,6 +41,11 @@ // column-major matrix, the solver will actually solve A^T x = b. // +extern Eigen::MatrixXd init_vertices; +extern Eigen::MatrixXd test_vertices; +extern std::vector test_boundary_nodes; +Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes); + namespace polysolve::linear { From 03fc84f5a6f1bc820449a548b028d56dc5c5d97b Mon Sep 17 00:00:00 2001 From: xDarkLemon Date: Tue, 19 Mar 2024 01:58:00 -0400 Subject: [PATCH 3/5] Fix compile for trilinos --- .github/workflows/continuous.yml | 5 ++-- CMakeLists.txt | 10 +++++--- cmake/recipes/trilinos.cmake | 42 ++++++++++++++++++-------------- 3 files changed, 34 insertions(+), 23 deletions(-) diff --git a/.github/workflows/continuous.yml b/.github/workflows/continuous.yml index 352d4ffb..2fc06bbb 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 c5d8cb08..aad24995 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,7 +77,7 @@ 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) @@ -288,8 +288,12 @@ 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) + 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 diff --git a/cmake/recipes/trilinos.cmake b/cmake/recipes/trilinos.cmake index 1b15890e..15174eae 100644 --- a/cmake/recipes/trilinos.cmake +++ b/cmake/recipes/trilinos.cmake @@ -7,31 +7,37 @@ 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() -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 ) +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 From 8ab068d656b6ea5e1b591508f1fd57f836f34ac7 Mon Sep 17 00:00:00 2001 From: xDarkLemon Date: Mon, 12 Aug 2024 16:29:34 -0400 Subject: [PATCH 4/5] Add Trilinos nullspace vector --- src/polysolve/linear/Solver.hpp | 2 + src/polysolve/linear/TrilinosSolver.cpp | 130 +++++++++++++++--------- src/polysolve/linear/TrilinosSolver.hpp | 6 +- 3 files changed, 87 insertions(+), 51 deletions(-) diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index 8b8702c7..9a53121b 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -38,6 +38,7 @@ namespace polysolve::linear public: // Shortcut alias typedef Eigen::VectorXd VectorXd; + typedef Eigen::MatrixXd MatrixXd; template using Ref = Eigen::Ref; @@ -127,6 +128,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 index 4be30094..1edd18ff 100644 --- a/src/polysolve/linear/TrilinosSolver.cpp +++ b/src/polysolve/linear/TrilinosSolver.cpp @@ -7,34 +7,6 @@ #include //////////////////////////////////////////////////////////////////////////////// -Eigen::MatrixXd test_vertices; -Eigen::MatrixXd init_vertices; -std::vector test_boundary_nodes; -Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes) -{ - // Remove boundary vertices - if (boundary_nodes.empty()) - { - return vertices; - } - else - { - std::vector order_nodes = boundary_nodes; - std::sort(order_nodes.begin(), order_nodes.end()); - Eigen::MatrixXd out_vertices; - std::vector 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 { TrilinosSolver::TrilinosSolver() @@ -159,35 +131,101 @@ namespace polysolve::linear } - void TrilinosSolver::solve(const Eigen::Ref rhs, Eigen::Ref result) + 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 (is_nullspace_) + if (nullspace.cols()==3) { - if (test_vertices.cols()==3) - { - reduced_vertices=remove_boundary_vertices(test_vertices,test_boundary_nodes); - MLList.set("null space: type","elasticity from coordinates"); - MLList.set("x-coordinates", reduced_vertices.col(0).data()); - MLList.set("y-coordinates", reduced_vertices.col(1).data()); - MLList.set("z-coordinates", reduced_vertices.col(2).data()); - MLList.set("aggregation: threshold",0.00); + 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); } - if (test_vertices.cols()==2) - { - reduced_vertices=remove_boundary_vertices(test_vertices,test_boundary_nodes); - MLList.set("null space: type","elasticity from coordinates"); - MLList.set("x-coordinates", reduced_vertices.col(0).data()); - MLList.set("y-coordinates", reduced_vertices.col(1).data()); - MLList.set("aggregation: threshold",0.00); - } + 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()); diff --git a/src/polysolve/linear/TrilinosSolver.hpp b/src/polysolve/linear/TrilinosSolver.hpp index 56d839c6..f1cf8975 100644 --- a/src/polysolve/linear/TrilinosSolver.hpp +++ b/src/polysolve/linear/TrilinosSolver.hpp @@ -41,11 +41,6 @@ // column-major matrix, the solver will actually solve A^T x = b. // -extern Eigen::MatrixXd init_vertices; -extern Eigen::MatrixXd test_vertices; -extern std::vector test_boundary_nodes; -Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector &boundary_nodes); - namespace polysolve::linear { @@ -78,6 +73,7 @@ namespace polysolve::linear // 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"; } From f920ea44a396d6467685d20d69f8ddba0a4af24b Mon Sep 17 00:00:00 2001 From: xDarkLemon Date: Mon, 12 Aug 2024 17:27:18 -0400 Subject: [PATCH 5/5] Fix Trilinos nullspace --- src/polysolve/linear/Solver.hpp | 3 --- src/polysolve/linear/TrilinosSolver.cpp | 10 ---------- src/polysolve/linear/TrilinosSolver.hpp | 2 -- 3 files changed, 15 deletions(-) diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index 9a53121b..3df9128f 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -4,9 +4,6 @@ #include -#include -#include -#include #define POLYSOLVE_DELETE_MOVE_COPY(Base) \ Base(Base &&) = delete; \ diff --git a/src/polysolve/linear/TrilinosSolver.cpp b/src/polysolve/linear/TrilinosSolver.cpp index 1edd18ff..4e943068 100644 --- a/src/polysolve/linear/TrilinosSolver.cpp +++ b/src/polysolve/linear/TrilinosSolver.cpp @@ -12,7 +12,6 @@ namespace polysolve::linear TrilinosSolver::TrilinosSolver() { precond_num_ = 0; -// #ifdef HAVE_MPI int done_already; MPI_Initialized(&done_already); if (!done_already) @@ -25,9 +24,6 @@ namespace polysolve::linear MPI_Init(&argc, &argvv); CommPtr = new Epetra_MpiComm(MPI_COMM_WORLD); } -// #else - // CommPtr=new Epetra_SerialComm; -// #endif } //////////////////////////////////////////////////////////////// @@ -50,10 +46,6 @@ namespace polysolve::linear { conv_tol_ = params["Trilinos"]["tolerance"]; } - if (params["Trilinos"].contains("is_nullspace")) - { - is_nullspace_ = params["Trilinos"]["is_nullspace"]; - } } } @@ -289,9 +281,7 @@ namespace polysolve::linear delete A; delete rowMap; delete MLPrec; -// #ifdef HAVE_MPI MPI_Finalize() ; -// #endif } } diff --git a/src/polysolve/linear/TrilinosSolver.hpp b/src/polysolve/linear/TrilinosSolver.hpp index f1cf8975..0b3fdb20 100644 --- a/src/polysolve/linear/TrilinosSolver.hpp +++ b/src/polysolve/linear/TrilinosSolver.hpp @@ -84,8 +84,6 @@ namespace polysolve::linear double conv_tol_ = 1e-8; size_t iterations_; double residual_error_; - bool is_nullspace_ = true; - Eigen::MatrixXd reduced_vertices; ML_Epetra::MultiLevelPreconditioner* MLPrec=NULL; Epetra_Map *rowMap=NULL;