branch: remove-local-superlu commit bb9888000da89f9d55f2f474545994e76f4063ae Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Sat Nov 18 20:52:37 2023 +0100
Make SuperLU optional --- interface/src/getfemint_precond.h | 43 ++++++++++++++++++++++++--------------- interface/src/gf_linsolve.cc | 5 ++++- interface/src/gf_precond.cc | 4 ++++ src/Makefile.am | 7 +++++-- src/getfem/getfem_config.h | 7 +++++++ src/getfem/getfem_model_solvers.h | 18 ++++++++++++---- src/getfem/getfem_superlu.h | 10 ++++----- tests/laplacian.cc | 4 +++- tests/schwarz_additive.cc | 28 +++++++++++++++---------- 9 files changed, 86 insertions(+), 40 deletions(-) diff --git a/interface/src/getfemint_precond.h b/interface/src/getfemint_precond.h index 8af6c79b..7dee1aaa 100644 --- a/interface/src/getfemint_precond.h +++ b/interface/src/getfemint_precond.h @@ -52,7 +52,7 @@ namespace getfemint { size_type ncols(void) const { return gsp ? gsp->ncols() : ncols_; } void set_dimensions(size_type m, size_type n) { nrows_ = m; ncols_ = n; } gprecond_base() : nrows_(0), ncols_(0), type(IDENTITY), gsp(0) {} - const char *name() const { + const char *name() const { const char *p[] = { "IDENTITY", "DIAG", "ILDLT", "ILDLTT", "ILU", "ILUT", "SUPERLU", "GSPARSE" }; return p[type]; @@ -69,7 +69,9 @@ namespace getfemint { std::unique_ptr<gmm::ildltt_precond<cscmat> > ildltt; std::unique_ptr<gmm::ilu_precond<cscmat> > ilu; std::unique_ptr<gmm::ilut_precond<cscmat> > ilut; +#if defined(GETFEM_USES_SUPERLU) std::unique_ptr<gmm::SuperLU_factor<T> > superlu; +#endif virtual size_type memsize() const { size_type sz = sizeof(*this); @@ -81,10 +83,15 @@ namespace getfemint { case ILDLT: sz += ildlt->memsize(); break; case ILDLTT: sz += ildltt->memsize(); break; case SUPERLU: - sz += size_type(superlu->memsize()); break; +#if defined(GETFEM_USES_SUPERLU) + sz += size_type(superlu->memsize()); +#else + GMM_ASSERT1(false, "GetFEM built without SuperLU support"); +#endif + break; case SPMAT: sz += gsp->memsize(); break; } - return sz; + return sz; } }; @@ -120,28 +127,32 @@ namespace gmm { switch (precond.type) { case getfemint::gprecond_base::IDENTITY: gmm::copy(v,w); break; case getfemint::gprecond_base::DIAG: - gmm::mult(*precond.diagonal, v, w); - break; - case getfemint::gprecond_base::ILDLT: - if (do_mult) gmm::mult(*precond.ildlt, v, w); + gmm::mult(*precond.diagonal, v, w); + break; + case getfemint::gprecond_base::ILDLT: + if (do_mult) gmm::mult(*precond.ildlt, v, w); else gmm::transposed_mult(*precond.ildlt, v, w); break; - case getfemint::gprecond_base::ILDLTT: - if (do_mult) gmm::mult(*precond.ildltt, v, w); + case getfemint::gprecond_base::ILDLTT: + if (do_mult) gmm::mult(*precond.ildltt, v, w); else gmm::transposed_mult(*precond.ildltt, v, w); break; - case getfemint::gprecond_base::ILU: - if (do_mult) gmm::mult(*precond.ilu, v, w); + case getfemint::gprecond_base::ILU: + if (do_mult) gmm::mult(*precond.ilu, v, w); else gmm::transposed_mult(*precond.ilu, v, w); break; - case getfemint::gprecond_base::ILUT: - if (do_mult) gmm::mult(*precond.ilut, v, w); + case getfemint::gprecond_base::ILUT: + if (do_mult) gmm::mult(*precond.ilut, v, w); else gmm::transposed_mult(*precond.ilut, v, w); break; case getfemint::gprecond_base::SUPERLU: - if (do_mult) precond.superlu->solve(w,v); - else precond.superlu->solve(w,v,gmm::SuperLU_factor<T>::LU_TRANSP); - break; +#if defined(GETFEM_USES_SUPERLU) + if (do_mult) precond.superlu->solve(w,v); + else precond.superlu->solve(w,v,gmm::SuperLU_factor<T>::LU_TRANSP); +#else + GMM_ASSERT1(false, "GetFEM built without SuperLU support"); +#endif + break; case getfemint::gprecond_base::SPMAT: precond.gsp->mult_or_transposed_mult(v, w, !do_mult); break; diff --git a/interface/src/gf_linsolve.cc b/interface/src/gf_linsolve.cc index a04adbe5..1b8f42ff 100644 --- a/interface/src/gf_linsolve.cc +++ b/interface/src/gf_linsolve.cc @@ -85,6 +85,7 @@ void iterative_gmm_solver(iterative_gmm_solver_type stype, else iterative_gmm_solver(stype, gsp, in, out, scalar_type()); } +#if defined(GETFEM_USES_SUPERLU) template <typename T> static void superlu_solver(gsparse &gsp, getfemint::mexargs_in& in, getfemint::mexargs_out& out, T) { @@ -96,6 +97,7 @@ superlu_solver(gsparse &gsp, if (out.remaining()) out.pop().from_scalar(rcond ? 1./rcond : 0.); } +#endif #if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H) template <typename T> static void @@ -178,6 +180,7 @@ void gf_linsolve(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { ); +#if defined(GETFEM_USES_SUPERLU) /*@FUNC @CELL{U, cond} = ('lu', @tsp M, @vec b) Alias for ::LINSOLVE('superlu',...)@*/ sub_command @@ -190,7 +193,6 @@ void gf_linsolve(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { else superlu_solver(gsp, in, out, scalar_type()); ); - /*@FUNC @CELL{U, cond} = ('superlu', @tsp M, @vec b) Solve `M.U = b` apply the SuperLU solver (sparse LU factorization). @@ -204,6 +206,7 @@ void gf_linsolve(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { if (gsp.is_complex()) superlu_solver(gsp, in, out, complex_type()); else superlu_solver(gsp, in, out, scalar_type()); ); +#endif #if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H) /*@FUNC @CELL{U, cond} = ('mumps', @tsp M, @vec b) diff --git a/interface/src/gf_precond.cc b/interface/src/gf_precond.cc index 970eeff7..7e6b0288 100644 --- a/interface/src/gf_precond.cc +++ b/interface/src/gf_precond.cc @@ -68,6 +68,7 @@ precond_ilut(gsparse &M, int additional_fillin, double threshold, mexargs_out& o p.ilut = std::make_unique<gmm::ilut_precond<typename gprecond<T>::cscmat>>(M.csc(T()), additional_fillin, threshold); } +#if defined(GETFEM_USES_SUPERLU) template <typename T> static void precond_superlu(gsparse &M, mexargs_out& out, T) { gprecond<T> &p = precond_new(out, T()); @@ -75,6 +76,7 @@ precond_superlu(gsparse &M, mexargs_out& out, T) { p.superlu = std::make_unique<gmm::SuperLU_factor<T>>(); p.superlu.get()->build_with(M.csc(T())); } +#endif static void precond_spmat(gsparse *gsp, mexargs_out& out) { if (gsp->is_complex()) { @@ -213,6 +215,7 @@ void gf_precond(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { out, scalar_type()); ); +#if defined(GETFEM_USES_SUPERLU) /*@INIT PC = ('superlu', @tsp m) Uses SuperLU to build an exact factorization of the sparse matrix `m`. This preconditioner is only available if the getfem-interface was @@ -224,6 +227,7 @@ void gf_precond(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { if (M->is_complex()) precond_superlu(*M, out, complex_type()); else precond_superlu(*M, out, scalar_type()); ); +#endif /*@INIT PC = ('spmat', @tsp m) Preconditioner given explicitely by a sparse matrix.@*/ diff --git a/src/Makefile.am b/src/Makefile.am index 16d06fcc..ee1268a7 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -158,7 +158,6 @@ nobase_include_HEADERS = \ getfem/getfem_nonlinear_elasticity.h \ getfem/getfem_fourth_order.h \ getfem/getfem_Navier_Stokes.h \ - getfem/getfem_superlu.h \ getfem/getfem_plasticity.h \ getfem/getfem_omp.h \ getfem/getfem_continuation.h \ @@ -190,7 +189,6 @@ SRC = \ bgeot_ftool.cc \ getfem_models.cc \ getfem_model_solvers.cc \ - getfem_superlu.cc \ getfem_mesh.cc \ getfem_mesh_region.cc \ getfem_context.cc \ @@ -247,6 +245,11 @@ SRC = \ getfem_continuation.cc # getfem_enumeration_dof_para.cc +if SUPERLU +nobase_include_HEADERS += getfem/getfem_superlu.h +SRC += getfem_superlu.cc +endif + getfem/getfem_im_list.h : ../cubature/getfem_im_list.h cp ../cubature/getfem_im_list.h getfem/getfem_im_list.h diff --git a/src/getfem/getfem_config.h b/src/getfem/getfem_config.h index a06821bc..27f4ad17 100644 --- a/src/getfem/getfem_config.h +++ b/src/getfem/getfem_config.h @@ -167,6 +167,10 @@ #define GETFEM_MPI_INIT(argc, argv) {GMM_TRACE1("Running sequential Getfem");} #define GETFEM_MPI_FINALIZE {} +#if defined(HAVE_SUPERLU_SLU_CNAMES_H) +# define GETFEM_USES_SUPERLU +#endif + #if defined(GETFEM_HAVE_DMUMPS_C_H) # ifndef GMM_USES_MUMPS # define GMM_USES_MUMPS @@ -211,7 +215,10 @@ #include "bgeot_tensor.h" #include "bgeot_poly.h" + +#if defined(GETFEM_USES_SUPERLU) #include "getfem_superlu.h" +#endif /// GEneric Tool for Finite Element Methods. namespace getfem { diff --git a/src/getfem/getfem_model_solvers.h b/src/getfem/getfem_model_solvers.h index 469896d8..e3397197 100644 --- a/src/getfem/getfem_model_solvers.h +++ b/src/getfem/getfem_model_solvers.h @@ -141,6 +141,7 @@ namespace getfem { } }; +#ifdef GETFEM_USES_SUPERLU template <typename MAT, typename VECT> struct linear_solver_superlu : public abstract_linear_solver<MAT, VECT> { @@ -155,6 +156,7 @@ namespace getfem { if (iter.get_noisy()) cout << "condition number: " << 1.0/rcond<< endl; } }; +#endif template <typename MAT, typename VECT> struct linear_solver_dense_lu : public abstract_linear_solver<MAT, VECT> { @@ -631,17 +633,20 @@ namespace getfem { <linear_solver_distributed_mumps<MATRIX, VECTOR>>(); #else size_type ndof = md.nb_dof(), max3d = 15000, dim = md.leading_dimension(); -# ifdef GMM_USES_MUMPS +# if defined(GMM_USES_MUMPS) max3d = 250000; # endif if ((ndof<300000 && dim<=2) || (ndof<max3d && dim<=3) || (ndof<1000)) { -# ifdef GMM_USES_MUMPS +# if defined(GMM_USES_MUMPS) if (md.is_symmetric()) return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>(); else return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>(); -# else +# elif definded(GETFEM_USES_SUPERLU) return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>(); +# else + static_assert(false, + "At least one direct solver (MUMPS or SuperLU) is required"); # endif } else { @@ -665,8 +670,13 @@ namespace getfem { std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>> select_linear_solver(const model &md, const std::string &name) { std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>> p; - if (bgeot::casecmp(name, "superlu") == 0) + if (bgeot::casecmp(name, "superlu") == 0) { +#ifdef GETFEM_USES_SUPERLU return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>(); +#else + GMM_ASSERT1(false, "SuperLU is not interfaced"); +#endif + } else if (bgeot::casecmp(name, "dense_lu") == 0) return std::make_shared<linear_solver_dense_lu<MATRIX, VECTOR>>(); else if (bgeot::casecmp(name, "mumps") == 0) { diff --git a/src/getfem/getfem_superlu.h b/src/getfem/getfem_superlu.h index 06b7c738..8e3255ce 100644 --- a/src/getfem/getfem_superlu.h +++ b/src/getfem/getfem_superlu.h @@ -38,12 +38,11 @@ does not include any of the superlu headers, hence when getfem is installed, it does not need to install the superlu headers. */ +#ifdef GETFEM_USES_SUPERLU + +#ifndef GETFEM_SUPERLU_H +#define GETFEM_SUPERLU_H -#ifndef GETFEM_SUPERLU -#define GETFEM_SUPERLU -#ifndef GMM_USES_SUPERLU -#define GMM_USES_SUPERLU -#endif #include "getfem_config.h" #include "gmm/gmm_kernel.h" @@ -128,3 +127,4 @@ namespace gmm { extern "C" void set_superlu_callback(int (*cb)()); #endif +#endif diff --git a/tests/laplacian.cc b/tests/laplacian.cc index d4caf948..370f5e7f 100644 --- a/tests/laplacian.cc +++ b/tests/laplacian.cc @@ -283,8 +283,10 @@ bool laplacian_problem::solve(void) { // gmm::cg(SM, U, B, P, iter); gmm::gmres(SM, U, B, P, 50, iter); } else { - double rcond; + double rcond; +#ifdef GETFEM_USES_SUPERLU gmm::SuperLU_solve(SM, U, B, rcond); +#endif cout << "cond = " << 1/rcond << "\n"; } diff --git a/tests/schwarz_additive.cc b/tests/schwarz_additive.cc index 9b58bd4e..abc957d3 100644 --- a/tests/schwarz_additive.cc +++ b/tests/schwarz_additive.cc @@ -26,8 +26,6 @@ /* */ /**************************************************************************/ -#define GMM_USES_SUPERLU - #include "getfem/getfem_assembling.h" #include "getfem/getfem_regular_meshes.h" #include "getfem/getfem_export.h" @@ -72,20 +70,24 @@ struct pb_data { linalg_vector U, F; /* Unknown and right hand side. */ int solver; - void assemble(void); + void assemble(); void init(bgeot::md_param ¶ms); - int solve_cg(void); - int solve_cg2(void); - int solve_superlu(void); + int solve_cg(); + int solve_cg2(); +#if defined(GETFEM_USES_SUPERLU) + int solve_superlu(); +#endif int solve_schwarz(int); - int solve(void) { + int solve() { cout << "solving" << endl; switch (solver) { case 0 : return solve_cg(); case 1 : return solve_cg2(); +#if defined(GETFEM_USES_SUPERLU) case 2 : return solve_superlu(); +#endif default : return solve_schwarz(solver); } return 0; @@ -190,7 +192,7 @@ void pb_data::init(bgeot::md_param ¶ms) { } } -void pb_data::assemble(void) { +void pb_data::assemble() { size_type nb_dof = mef.nb_dof(); std::cout << "number of dof : "<< nb_dof << endl; size_type nb_dof_data = mef_data.nb_dof(); @@ -217,20 +219,22 @@ void pb_data::assemble(void) { getfem::assembling_Dirichlet_condition(RM, F, mef, 0, UD); } -int pb_data::solve_cg(void) { +int pb_data::solve_cg() { gmm::iteration iter(residual, 1, 1000000); gmm::ildlt_precond<general_sparse_matrix> P(RM); gmm::cg(RM, U, F, gmm::identity_matrix(), P, iter); return int(iter.get_iteration()); } -int pb_data::solve_superlu(void) { +#if defined(GETFEM_USES_SUPERLU) +int pb_data::solve_superlu() { double rcond; SuperLU_solve(RM, U, F, rcond); return 1; } +#endif -int pb_data::solve_cg2(void) { +int pb_data::solve_cg2() { gmm::iteration iter(residual, 1, 1000000); gmm::cg(RM, U, F, gmm::identity_matrix(), gmm::identity_matrix(), iter); return int(iter.get_iteration()); @@ -269,10 +273,12 @@ int pb_data::solve_schwarz(int version) { gmm::ilu_precond<general_sparse_matrix>(), vB, iter, gmm::using_gmres(), gmm::using_gmres()); break; +#if defined(GETFEM_USES_SUPERLU) case 5 : gmm::additive_schwarz(RM, U, F, gmm::ilu_precond<general_sparse_matrix>(), vB, iter, gmm::using_superlu(), gmm::using_cg()); break; +#endif } return 0; }