[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: |
Tue, 23 Jan 2024 05:31:40 -0500 (EST) |
branch: add-umfpack-interface
commit a3afc92c2b6a10b138da320b28e06baa75746e94
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Tue Jan 23 11:30:21 2024 +0100
Add experimental UMFPACK support
- compiles but does not work properly yet
---
CMakeLists.txt | 1 +
cmake/gmm_arch_config.h.in | 3 +
configure.ac | 107 ++++++++++++++++++-
src/Makefile.am | 1 +
src/getfem/getfem_model_solvers.h | 21 ++++
src/gmm/gmm_UMFPACK_interface.h | 188 ++++++++++++++++++++++++++++++++++
src/gmm/gmm_arch_config.h.in | 3 +
src/gmm/gmm_solver_Schwarz_additive.h | 8 +-
tests/test_condensation.cc | 4 +-
9 files changed, 329 insertions(+), 7 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7ff2cdd8..be062820 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -322,6 +322,7 @@ set(HEADERS
src/gmm/gmm_superlu_interface.h
src/gmm/gmm_transposed.h
src/gmm/gmm_tri_solve.h
+ src/gmm/gmm_UMFPACK_interface.h
src/gmm/gmm_vector.h
src/gmm/gmm_vector_to_matrix.h
src/getfem/bgeot_comma_init.h
diff --git a/cmake/gmm_arch_config.h.in b/cmake/gmm_arch_config.h.in
index 50f47daa..4b4b029a 100644
--- a/cmake/gmm_arch_config.h.in
+++ b/cmake/gmm_arch_config.h.in
@@ -19,6 +19,9 @@
/* defined if GMM is linked to the superlu library */
#cmakedefine GMM_USES_SUPERLU
+/* defined if GMM is linked to the umfpack library */
+#cmakedefine GMM_USES_UMFPACK
+
/* Use blas with 64 bits integers */
#cmakedefine GMM_USE_BLAS64_INTERFACE
diff --git a/configure.ac b/configure.ac
index cd519ed5..65e505ce 100644
--- a/configure.ac
+++ b/configure.ac
@@ -958,9 +958,98 @@ fi
dnl ---------------------------END OF MUMPS TEST--------------------------
+dnl ------------------------------UMFPACK TEST----------------------------
+require_umfpack="auto"
+AC_ARG_ENABLE(umfpack,
+ [AS_HELP_STRING([--enable-umfpack], [Enable UMFPACK support])],
+ [require_umfpack=$enableval],
+ [require_umfpack="auto"])
+
+UMFPACK_LIBS="-lumfpack"
+# the user can override these defaults using --with-umfpack=
+AC_ARG_WITH(umfpack,
+ [AS_HELP_STRING([--with-umfpack=<lib>],[use UMFPACK library <lib>])],
+ [case $with_umfpack in
+ yes | "")
+ if test "x$require_umfpack" = "xno"; then
+ AC_MSG_ERROR([Contradicting arguments between --enable-umfpack and
--with-umfpack.])
+ elif test "x$require_umfpack" = "xauto"; then
+ require_umfpack="yes"
+ fi;;
+ no)
+ if test "x$require_umfpack" = "xyes"; then
+ AC_MSG_ERROR([Contradicting arguments between --enable-umfpack and
--with-umfpack.])
+ elif test "x$require_umfpack" = "xauto"; then
+ require_umfpack="no"
+ fi;;
+ -* | */* | *.a | *.so | *.so.* | *.o| builtin)
UMFPACK_LIBS="$with_umfpack";;
+ *) UMFPACK_LIBS=`echo $with_umfpack | sed -e 's/^/-l/g;s/ \+/ -l/g'`;;
+ esac]
+)
+
+UMFPACKINC=""
+AC_ARG_WITH(umfpack-include-dir,
+ [AS_HELP_STRING([--with-umfpack-include-dir],[directory in which the
umfpack.h headers can be found])],
+ [ if test x$require_umfpack = xno; then
+ AC_MSG_ERROR([Inconsistent options for --enable-umfpack, --with-umfpack
and --with-umfpack-include-dir.]);
+ else
+ require_umfpack="yes"
+ case $withval in
+ -I* ) UMFPACKINC="$withval";;
+ * ) UMFPACKINC="-I$withval";;
+ esac
+ fi;],
+)
+CPPFLAGS="$CPPFLAGS $UMFPACKINC"
-if test "x$found_superlu" = "xno" -a "x$found_mumps" = "xno"; then
- AC_MSG_ERROR([Neither MUMPS nor SuperLU was enabled. At least one linear
solver is required.])
+if test "x$require_umfpack" = "xno"; then
+ UMFPACK_LIBS=""
+ found_umfpack="no"
+ echo "Building with UMFPACK explicitly disabled";
+else
+ AC_CHECK_HEADERS(
+ [umfpack.h],
+ [found_umfpack="yes"],
+ [ if test "x$require_umfpack" = "xyes"; then
+ AC_MSG_ERROR([Header files of UMFPACK not found.]);
+ else
+ found_umfpack="no"
+ fi;
+ ])
+ if test x$found_umfpack = xyes; then
+ save_LIBS="$LIBS";
+ AC_CHECK_LIB([umfpack], [umfpack_di_solve],
+ [AC_DEFINE(GMM_USES_UMFPACK,,[defined if GMM is linked to the
umfpack library])],
+ [if test "x$require_umfpack" = "xyes"; then
+ AC_MSG_ERROR([UMFPACK library not found]);
+ else
+ found_umfpack="no"
+ fi;])
+ if test "x$found_umfpack" = "xyes"; then
+ echo "Building with UMFPACK (use --enable-umfpack=no to disable it)"
+ LIBS="$UMFPACK_LIBS $save_LIBS"
+ else
+ UMFPACK_LIBS=""
+ LIBS="$save_LIBS"
+ fi
+ elif test "x$require_umfpack" = "xyes"; then
+ AC_MSG_ERROR([UMFPACK header files not found but required by the user.
Aborting configure...]);
+ else
+ echo "UMFPACK header files not found, building without UMFPACK"
+ UMPFPACK_LIBS=""
+ fi
+fi
+
+AM_CONDITIONAL(UMFPACK, test x$found_umfpack = xyes)
+AC_SUBST([UMFPACK_LIBS])
+if test "x$found_umfpack" = "xyes"; then
+ echo "Configuration of UMFPACK done"
+fi
+dnl ---------------------------END OF UMFPACK TEST------------------------
+
+
+if test "x$found_superlu" = "xno" -a "x$found_mumps" = "xno" -a
"x$found_umfpack" = "xno"; then
+ AC_MSG_ERROR([Neither MUMPS, SuperLU, nor UMFPACK was enabled. At least one
external direct linear solver is required.])
exit 1
fi
@@ -1273,12 +1362,22 @@ else
fi;
if test "x$found_mumps" = "xyes"; then
- echo "- Mumps found. A direct solver for large sparse linear systems."
+ echo "- MUMPS found. A direct solver for large sparse linear systems."
else
if test "x$require_mumps" = "xno"; then
echo "- Not using the MUMPS library for large sparse linear systems."
else
- echo "- Mumps not found. Not using the MUMPS library for large sparse
linear systems."
+ echo "- MUMPS not found. Not using the MUMPS library for large sparse
linear systems."
+ fi
+fi;
+
+if test "x$found_umfpack" = "xyes"; then
+ echo "- UMFPACK found. A direct solver for large sparse linear systems."
+else
+ if test "x$require_umfpack" = "xno"; then
+ echo "- Not using the UMFPACK library for large sparse linear systems."
+ else
+ echo "- UMFPACK not found. Not using the UMFPACK library for large sparse
linear systems."
fi
fi;
diff --git a/src/Makefile.am b/src/Makefile.am
index c002f042..61b3da0c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -73,6 +73,7 @@ nobase_include_HEADERS = \
gmm/gmm_except.h \
gmm/gmm_feedback_management.h \
gmm/gmm_MUMPS_interface.h \
+ gmm/gmm_UMFPACK_interface.h \
getfem/dal_config.h \
getfem/dal_singleton.h \
getfem/dal_basic.h \
diff --git a/src/getfem/getfem_model_solvers.h
b/src/getfem/getfem_model_solvers.h
index da570c01..01112a6c 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -42,6 +42,7 @@
#include "getfem_models.h"
#include "gmm/gmm_superlu_interface.h"
#include "gmm/gmm_MUMPS_interface.h"
+#include "gmm/gmm_UMFPACK_interface.h"
#include "gmm/gmm_iter.h"
#include "gmm/gmm_iter_solvers.h"
#include "gmm/gmm_dense_qr.h"
@@ -190,6 +191,17 @@ namespace getfem {
};
#endif
+#if defined(GMM_USES_UMFPACK)
+ template <typename MAT, typename VECT>
+ struct linear_solver_umfpack : public abstract_linear_solver<MAT, VECT> {
+ void operator ()(const MAT &M, VECT &x, const VECT &b,
+ gmm::iteration &iter) const {
+ bool ok = gmm::UMFPACK_solve(M, x, b);
+ iter.enforce_converged(ok);
+ }
+ };
+#endif
+
#if GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
template <typename MAT, typename VECT>
struct linear_solver_distributed_mumps
@@ -645,6 +657,8 @@ namespace getfem {
return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
# elif defined(GMM_USES_SUPERLU)
return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>();
+# elif defined(GMM_USES_UMFPACK)
+ return std::make_shared<linear_solver_umfpack<MATRIX, VECTOR>>();
# else
static_assert(false,
"At least one direct solver (MUMPS or SuperLU) is
required");
@@ -690,6 +704,13 @@ namespace getfem {
# endif
#else
GMM_ASSERT1(false, "Mumps is not interfaced");
+#endif
+ }
+ else if (bgeot::casecmp(name, "umfpack") == 0) {
+#if defined(GMM_USES_UMFPACK)
+ return std::make_shared<linear_solver_umfpack<MATRIX, VECTOR>>();
+#else
+ GMM_ASSERT1(false, "UMFPACK is not interfaced");
#endif
}
else if (bgeot::casecmp(name, "cg/ildlt") == 0)
diff --git a/src/gmm/gmm_UMFPACK_interface.h b/src/gmm/gmm_UMFPACK_interface.h
new file mode 100644
index 00000000..1efe2602
--- /dev/null
+++ b/src/gmm/gmm_UMFPACK_interface.h
@@ -0,0 +1,188 @@
+/* -*- c++ -*- (enables emacs c++ mode) */
+/*===========================================================================
+
+ Copyright (C) 2024-2024 Konstantinos Poulios
+
+ This file is a part of GetFEM
+
+ GetFEM is free software; you can redistribute it and/or modify it
+ under the terms of the GNU Lesser General Public License as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version along with the GCC Runtime Library
+ Exception either version 3.1 or (at your option) any later version.
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+ License and GCC Runtime Library Exception for more details.
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+ As a special exception, you may use this file as it is a part of a free
+ software library without restriction. Specifically, if other files
+ instantiate templates or use macros or inline functions from this file,
+ or you compile this file and link it with other files to produce an
+ executable, this file does not by itself cause the resulting executable
+ to be covered by the GNU Lesser General Public License. This exception
+ does not however invalidate any other reasons why the executable file
+ might be covered by the GNU Lesser General Public License.
+
+===========================================================================*/
+
+/**@file gmm_UMFPACK_interface.h
+ @author Konstantinos Poulios <logari81@gmail.com>,
+ @date January 11, 2024.
+ @brief Interface with UMFPACK (direct solver for sparse matrices).
+*/
+#if defined(GMM_USES_UMFPACK)
+
+#ifndef GMM_UMFPACK_INTERFACE_H
+#define GMM_UMFPACK_INTERFACE_H
+
+#include "gmm_kernel.h"
+
+
+extern "C" {
+#include <umfpack.h>
+}
+
+namespace gmm {
+
+ /* ********************************************************************* */
+ /* UMFPACK solve interface */
+ /* ********************************************************************* */
+
+ //
+ template <typename T>
+ struct umfpack_interf {
+ int n;
+ std::vector<int> c;
+ std::vector<int> i;
+ std::vector<double> ar;
+ std::vector<double> xr;
+ std::vector<double> br;
+
+ inline int
+ umfpack_symbolic(void **output) {
+ return umfpack_di_symbolic(n, n, &(c[0]), &(i[0]), &(ar[0]),
+ output, NULL, NULL);
+ }
+ inline int
+ umfpack_numeric(void *symbolic, void **output) {
+ return umfpack_di_numeric(&(c[0]), &(i[0]), &(ar[0]),
+ symbolic, output, NULL, NULL);
+ }
+ static inline void
+ umfpack_free_symbolic(void **symbolic) {
+ umfpack_di_free_symbolic(symbolic);
+ }
+ inline int
+ umfpack_solve(void *numeric) {
+ return umfpack_di_solve(UMFPACK_A, &(c[0]), &(i[0]), &(ar[0]),
+ &(xr[0]), &(br[0]), numeric, NULL, NULL);
+ }
+ static inline void
+ umfpack_free_numeric(void **numeric) {
+ umfpack_di_free_numeric(numeric);
+ }
+ template<typename MAT, typename VECTX, typename VECTB>
+ umfpack_interf(MAT &A, VECTX &X, VECTB &B)
+ {
+ n = int(gmm::vect_size(B));
+ GMM_ASSERT2(gmm::mat_nrows(A) == size_type(n), "Incompatible
matrix-vector sizes");
+ GMM_ASSERT2(gmm::mat_ncols(A) == size_type(n), "System matrix needs to
be square");
+ xr.resize(n);
+ br.resize(n);
+ gmm::copy(X, xr);
+ gmm::copy(B, br);
+ csc_matrix<double> csc_A(n, n);
+ gmm::copy(A, csc_A);
+ c.resize(csc_A.jc.size());
+ i.resize(csc_A.ir.size());
+ ar.resize(csc_A.pr.size());
+ gmm::copy(csc_A.jc, c);
+ gmm::copy(csc_A.ir, i);
+ gmm::copy(csc_A.pr, ar);
+ }
+ };
+
+ template <typename T>
+ struct umfpack_interf<std::complex<T>> {
+ int n;
+ std::vector<int> c;
+ std::vector<int> i;
+ std::vector<double> ar, ai;
+ std::vector<double> xr, xi;
+ std::vector<double> br, bi;
+
+ inline int
+ umfpack_symbolic(void **output) {
+ return umfpack_zi_symbolic(n, n, &(c[0]), &(i[0]), &(ar[0]), &(ai[0]),
+ output, NULL, NULL);
+ }
+ inline int
+ umfpack_numeric(void *symbolic, void **output) {
+ return umfpack_zi_numeric(&(c[0]), &(i[0]), &(ar[0]), &(ai[0]),
+ symbolic, output, NULL, NULL);
+ }
+ static inline void
+ umfpack_free_symbolic(void **symbolic) {
+ umfpack_zi_free_symbolic(symbolic);
+ }
+ inline int
+ umfpack_solve(void *numeric) {
+ return umfpack_zi_solve(UMFPACK_A, &(c[0]), &(i[0]), &(ar[0]), &(ai[0]),
+ &(xr[0]), &(xi[0]), &(br[0]), &(bi[0]), numeric,
NULL, NULL);
+ }
+ static inline void
+ umfpack_free_numeric(void **numeric) {
+ umfpack_zi_free_numeric(numeric);
+ }
+ template<typename MAT, typename VECTX, typename VECTB>
+ umfpack_interf(MAT &A, VECTX &X, VECTB &B)
+ {
+ n = int(gmm::vect_size(B));
+ GMM_ASSERT2(gmm::mat_nrows(A) == size_type(n), "Incompatible
matrix-vector sizes");
+ GMM_ASSERT2(gmm::mat_ncols(A) == size_type(n), "System matrix needs to
be square");
+ xr.resize(n);
+ xi.resize(n);
+ br.resize(n);
+ bi.resize(n);
+ gmm::copy(gmm::real_part(X), xr);
+ gmm::copy(gmm::imag_part(X), xi);
+ gmm::copy(gmm::real_part(B), br);
+ gmm::copy(gmm::imag_part(B), bi);
+ csc_matrix<std::complex<double>> csc_A(n, n);
+ gmm::copy(A, csc_A);
+ c.resize(csc_A.jc.size());
+ i.resize(csc_A.ir.size());
+ ar.resize(csc_A.pr.size());
+ ai.resize(csc_A.pr.size());
+ gmm::copy(csc_A.jc, c);
+ gmm::copy(csc_A.ir, i);
+ gmm::copy(gmm::real_part(csc_A.pr), ar);
+ gmm::copy(gmm::imag_part(csc_A.pr), ai);
+ }
+ };
+
+ // UMFPACK solve interface, returns "true" if successful
+ template <typename MAT, typename VECTX, typename VECTB>
+ bool UMFPACK_solve(const MAT &A, VECTX &X, const VECTB &B) {
+ typedef typename linalg_traits<MAT>::value_type T;
+ umfpack_interf<T> intrf(A, B, X); // possibly convert float to double
+ // umfpack supports only double precision
+ int status;
+ void *symbolic, *numeric;
+ status = intrf.umfpack_symbolic(&symbolic);
+ status = intrf.umfpack_numeric(symbolic, &numeric);
+ intrf.umfpack_free_symbolic(&symbolic);
+ status = intrf.umfpack_solve(numeric);
+ intrf.umfpack_free_numeric(&numeric);
+ return (status == UMFPACK_OK);
+ }
+
+}
+
+#endif // GMM_UMFPACK_INTERFACE_H
+
+#endif // GMM_USES_UMFPACK
diff --git a/src/gmm/gmm_arch_config.h.in b/src/gmm/gmm_arch_config.h.in
index 0d46b9ab..2c7dc709 100644
--- a/src/gmm/gmm_arch_config.h.in
+++ b/src/gmm/gmm_arch_config.h.in
@@ -18,6 +18,9 @@
/* defined if GMM is linked to the superlu library */
#undef GMM_USES_SUPERLU
+/* defined if GMM is linked to the umfpack library */
+#undef GMM_USES_UMFPACK
+
/* Use blas with 64 bits integers */
#undef GMM_USE_BLAS64_INTERFACE
diff --git a/src/gmm/gmm_solver_Schwarz_additive.h
b/src/gmm/gmm_solver_Schwarz_additive.h
index 378534e1..5a2c9cbc 100644
--- a/src/gmm/gmm_solver_Schwarz_additive.h
+++ b/src/gmm/gmm_solver_Schwarz_additive.h
@@ -41,8 +41,10 @@
#include "gmm_kernel.h"
#if defined(GMM_USES_SUPERLU)
#include "gmm_superlu_interface.h"
-#else
+#elif defined(GMM_USES_MUMPS)
#include "gmm_MUMPS_interface.h"
+#else //defined(GMM_USES_UMFPACK)
+#include "gmm_UMFPACK_interface.h"
#endif
#include "gmm_solver_cg.h"
#include "gmm_solver_gmres.h"
@@ -575,8 +577,10 @@ namespace gmm {
#if defined(GMM_USES_SUPERLU)
double rcond;
gmm::SuperLU_solve(M.vMloc[i], x, w, rcond);
-#else
+#elif defined(GMM_USES_MUMPS)
gmm::MUMPS_solve(M.vMloc[i], x, w);
+#else // defined(GMM_USES_UMFPACK)
+ gmm::UMFPACK_solve(M.vMloc[i], x, w);
#endif
// gmm::iteration iter(1E-10, 0, 100000);
//gmm::gmres(M.vMloc[i], x, w, gmm::identity_matrix(), 50, iter);
diff --git a/tests/test_condensation.cc b/tests/test_condensation.cc
index d62c3ed2..9509f53d 100644
--- a/tests/test_condensation.cc
+++ b/tests/test_condensation.cc
@@ -35,8 +35,10 @@ int main(int argc, char *argv[]) {
#if defined(GMM_USES_MUMPS)
std::string lsolver("mumps");
-#else
+#elif defined(GMM_USES_SUPERLU)
std::string lsolver("superlu");
+#else
+ std::string lsolver("umfpack");
#endif
gmm::set_traces_level(1);