branch: add-umfpack-interface commit a3afc92c2b6a10b138da320b28e06baa75746e94 Author: Konstantinos Poulios <logar...@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 <logar...@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);