[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Mon, 18 Dec 2023 11:06:20 -0500 (EST) |
branch: remove-local-superlu
commit aa84a361909484226d0d1c59f382b8c88d8f4c34
Author: Konstantinos Poulios <logari81@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 2df0dac8..1a997645 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;
}