From 9329ac73f6e3d3e6f3cddc820370d764dcb54ab4 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 26 Sep 2017 20:40:48 +0200 Subject: [PATCH] add umfpack solver to SparseMatrix --- cmake/Modules/FindIFEMDeps.cmake | 12 ++++++- cmake/Modules/FindISTL.cmake | 7 +--- src/IFEM.C | 11 +++++++ src/LinAlg/SparseMatrix.C | 55 ++++++++++++++++++++++++++++---- src/LinAlg/SparseMatrix.h | 7 +++- src/LinAlg/SystemMatrix.C | 1 + src/LinAlg/SystemMatrix.h | 2 +- src/SIM/SIMoptions.C | 4 +++ 8 files changed, 84 insertions(+), 15 deletions(-) diff --git a/cmake/Modules/FindIFEMDeps.cmake b/cmake/Modules/FindIFEMDeps.cmake index 93b0c804..02b1b41c 100644 --- a/cmake/Modules/FindIFEMDeps.cmake +++ b/cmake/Modules/FindIFEMDeps.cmake @@ -212,7 +212,17 @@ IF(IFEM_USE_OPENMP) ENDIF(OPENMP_FOUND) ENDIF(IFEM_USE_OPENMP) -# ISTL +# UMFPack +if(IFEM_USE_UMFPACK) + find_package(SuiteSparse COMPONENTS umfpack) + if(SuiteSparse_UMFPACK_FOUND) + list(APPEND IFEM_DEFINITIONS -DHAS_UMFPACK=1) + list(APPEND IFEM_DEPINCLUDES ${UMFPACK_INCLUDE_DIR}) + list(APPEND IFEM_DEPLIBS ${UMFPACK_LIBRARY}) + endif() +endif() + +# ISTL - need to be after UMFPack and SuperLU if(IFEM_USE_ISTL) find_package(ISTL) if(ISTL_FOUND) diff --git a/cmake/Modules/FindISTL.cmake b/cmake/Modules/FindISTL.cmake index 2aee35aa..dd66f952 100644 --- a/cmake/Modules/FindISTL.cmake +++ b/cmake/Modules/FindISTL.cmake @@ -20,17 +20,12 @@ if(SUPERLU_FOUND) list(APPEND ISTL_DEFINITIONS -DHAVE_SUPERLU=1 -DSUPERLU_POST_2005_VERSION=1 -DSUPERLU_NTYPE=1 -DSUPERLU_MIN_VERSION_4_3=1) endif() -if(IFEM_USE_UMFPACK) - find_package(SuiteSparse COMPONENTS umfpack) -endif() - -if(SuiteSparse_UMFPACK_FOUND) +if(SuiteSparse_UMFPACK_FOUND AND ISTL_FOUND) list(APPEND ISTL_DEFINITIONS -DHAVE_UMFPACK=1) list(APPEND ISTL_INCLUDE_DIRS ${UMFPACK_INCLUDE_DIR}) list(APPEND ISTL_LIBRARIES ${UMFPACK_LIBRARY}) endif() - list(APPEND ISTL_INCLUDE_DIRS ${ISTL_INCLUDEDIR}) include(FindPackageHandleStandardArgs) diff --git a/src/IFEM.C b/src/IFEM.C index c16aad58..0ab0f78d 100644 --- a/src/IFEM.C +++ b/src/IFEM.C @@ -24,6 +24,9 @@ #elif defined(HAVE_MPI) #include #endif +#ifdef HAS_UMFPACK +#include +#endif int IFEM::argc; char** IFEM::argv; @@ -113,6 +116,14 @@ int IFEM::Init (int arg_c, char** arg_v, const char* title) #else std::cout <<"disabled"; #endif + std::cout <<"\n UMFPack support: "; +#if HAS_UMFPACK + std::cout <<"enabled (v" << UMFPACK_MAIN_VERSION <<"." + << UMFPACK_SUB_VERSION <<"." + << UMFPACK_SUBSUB_VERSION << ")"; +#else + std::cout <<"disabled"; +#endif std::cout <<"\n ISTL support: "; #if HAS_ISTL diff --git a/src/LinAlg/SparseMatrix.C b/src/LinAlg/SparseMatrix.C index 162af2c3..dfad6136 100644 --- a/src/LinAlg/SparseMatrix.C +++ b/src/LinAlg/SparseMatrix.C @@ -22,6 +22,9 @@ #ifdef HAS_SAMG #include "samg.h" #endif +#ifdef HAS_UMFPACK +#include +#endif #ifdef USE_OPENMP #include #endif @@ -252,7 +255,7 @@ Real& SparseMatrix::operator () (size_t r, size_t c) return value; } } - else if (solver == SUPERLU) { + else if (solver == SUPERLU || solver == UMFPACK) { // Column-oriented format with 0-based indices IntVec::const_iterator begin = JA.begin() + IA[c-1]; IntVec::const_iterator end = JA.begin() + IA[c]; @@ -289,7 +292,7 @@ const Real& SparseMatrix::operator () (size_t r, size_t c) const ValueIter vit = elem.find(IJPair(r,c)); if (vit != elem.end()) return vit->second; } - else if (solver == SUPERLU) { + else if (solver == SUPERLU || solver == UMFPACK) { // Column-oriented format with 0-based indices IntVec::const_iterator begin = JA.begin() + IA[c-1]; IntVec::const_iterator end = JA.begin() + IA[c]; @@ -321,7 +324,7 @@ void SparseMatrix::dump (std::ostream& os, char format, const char* label) for (ValueIter it = elem.begin(); it != elem.end(); it++) os << it->first.first <<' '<< it->first.second <<" "<< it->second <<";\n"; - else if (solver == SUPERLU) { + else if (solver == SUPERLU || solver == UMFPACK) { // Column-oriented format with 0-based indices os << JA[0]+1 <<" 1 "<< A[0]; for (size_t j = 1; j <= ncol; j++) @@ -383,7 +386,7 @@ void SparseMatrix::printSparsity (std::ostream& os) const for (c = 1; c <= ncol; c++) if (editable) os << (elem.find(IJPair(r,c)) == elem.end() ? '.' : 'X'); - else if (solver == SUPERLU) { + else if (solver == SUPERLU || solver == UMFPACK) { // Column-oriented format with 0-based indices IntVec::const_iterator begin = JA.begin() + IA[c-1]; IntVec::const_iterator end = JA.begin() + IA[c]; @@ -488,7 +491,7 @@ bool SparseMatrix::add (const SystemMatrix& B, Real alpha) } else if (editable == 'P') { - if (solver == SUPERLU) + if (solver == SUPERLU || solver == UMFPACK) // Column-oriented format with 0-based indices for (size_t j = 1; j <= Bptr->ncol; j++) for (int i = Bptr->IA[j-1]; i < Bptr->IA[j]; i++) @@ -704,6 +707,7 @@ void SparseMatrix::preAssemble (const SAM& sam, bool delayLocking) << nrow <<"x"<< ncol <<"): "<< std::flush; switch (solver) { + case UMFPACK: case SUPERLU: this->optimiseSLU(dofc); break; case S_A_M_G: this->optimiseSAMG(); break; default: break; @@ -733,6 +737,7 @@ void SparseMatrix::preAssemble (const std::vector& MMNPC, size_t nel) } switch (solver) { + case UMFPACK: case SUPERLU: this->optimiseSLU(); break; case S_A_M_G: this->optimiseSAMG(); break; default: break; @@ -962,6 +967,7 @@ bool SparseMatrix::solve (SystemVector& B, bool, Real* rc) { case SUPERLU: return this->solveSLUx(*Bptr,rc); case S_A_M_G: return this->solveSAMG(*Bptr); + case UMFPACK: return this->solveUMF(*Bptr,rc); default: std::cerr <<"SparseMatrix::solve: No equation solver"<< std::endl; } @@ -1219,6 +1225,43 @@ bool SparseMatrix::solveSLUx (Vector& B, Real* rcond) } +bool SparseMatrix::solveUMF (Vector& B, Real* rcond) +{ + if (!factored) this->optimiseSLU(); + +#ifdef HAS_UMFPACK + double info[UMFPACK_INFO]; + Vector X(B.size()); + void* symbolic; + umfpack_di_symbolic(nrow, ncol, IA.data(), JA.data(), + A.data(), &symbolic, nullptr, info); + if (info[UMFPACK_STATUS] != UMFPACK_OK) + return false; + void* numeric; + umfpack_di_numeric(IA.data(), JA.data(), A.data(), symbolic, + &numeric, nullptr, info); + if (info[UMFPACK_STATUS] != UMFPACK_OK) + { + umfpack_di_free_symbolic(&symbolic); + return false; + } + if (rcond) + *rcond = info[UMFPACK_RCOND]; + umfpack_di_solve(UMFPACK_A, + IA.data(), JA.data(), A.data(), + &X[0], &B[0], numeric, nullptr, info); + if (info[UMFPACK_STATUS] == UMFPACK_OK) + B = X; + umfpack_di_free_symbolic(&symbolic); + umfpack_di_free_numeric(&numeric); + return info[UMFPACK_STATUS] == UMFPACK_OK; +#else + std::cerr <<"SparseMatrix::solve: UMFPACK solver not available"<< std::endl; + return false; +#endif +} + + bool SparseMatrix::solveSAMG (Vector& B) { int ierr = 1; @@ -1290,7 +1333,7 @@ Real SparseMatrix::Linfnorm () const if (editable) for (ValueIter it = elem.begin(); it != elem.end(); ++it) sums[it->first.first-1] += fabs(it->second); - else if (solver == SUPERLU) + else if (solver == SUPERLU || solver == UMFPACK) // Column-oriented format with 0-based row-indices for (size_t j = 1; j <= ncol; j++) for (int i = IA[j-1]; i < IA[j]; i++) diff --git a/src/LinAlg/SparseMatrix.h b/src/LinAlg/SparseMatrix.h index a5a1be3c..2e17a00d 100644 --- a/src/LinAlg/SparseMatrix.h +++ b/src/LinAlg/SparseMatrix.h @@ -39,7 +39,7 @@ class SparseMatrix : public SystemMatrix { public: //! \brief Available equation solvers for this matrix type. - enum SparseSolver { NONE, SUPERLU, S_A_M_G }; + enum SparseSolver { NONE, SUPERLU, S_A_M_G, UMFPACK }; //! \brief Default constructor creating an empty matrix. SparseMatrix(SparseSolver eqSolver = NONE, int nt = 1); @@ -236,6 +236,11 @@ protected: //! \param[out] rcond Reciprocal condition number of the LHS-matrix (optional) bool solveSLUx(Vector& B, Real* rcond); + //! \brief Invokes the UMFPACK equation solver for a given right-hand-side. + //! \param B Right-hand-side vector on input, solution vector on output + //! \param[out] rcond Reciprocal condition number of the LHS-matrix (optional) + bool solveUMF(Vector& B, Real* rcond); + //! \brief Writes the system matrix to the given output stream. virtual std::ostream& write(std::ostream& os) const; diff --git a/src/LinAlg/SystemMatrix.C b/src/LinAlg/SystemMatrix.C index 8e94c8fa..df4bfa89 100644 --- a/src/LinAlg/SystemMatrix.C +++ b/src/LinAlg/SystemMatrix.C @@ -118,6 +118,7 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType, case SPR : return new SPRMatrix(); case SPARSE: return new SparseMatrix(SparseMatrix::SUPERLU,num_thread_SLU); case SAMG : return new SparseMatrix(SparseMatrix::S_A_M_G); + case UMFPACK: return new SparseMatrix(SparseMatrix::UMFPACK); #ifdef HAS_ISTL case ISTL : return new ISTLMatrix(padm,defaultPar,ltype); #endif diff --git a/src/LinAlg/SystemMatrix.h b/src/LinAlg/SystemMatrix.h index a04478ea..84dedb08 100644 --- a/src/LinAlg/SystemMatrix.h +++ b/src/LinAlg/SystemMatrix.h @@ -203,7 +203,7 @@ class SystemMatrix public: //! \brief The available system matrix formats. enum Type { DENSE = 0, SPR = 1, SPARSE = 2, SAMG = 3, - PETSC = 4, ISTL = 5 }; + PETSC = 4, ISTL = 5, UMFPACK = 6 }; //! \brief Static method creating a matrix of the given type. static SystemMatrix* create(const ProcessAdm& padm, Type matrixType, diff --git a/src/SIM/SIMoptions.C b/src/SIM/SIMoptions.C index 5f869e4a..918d99a7 100644 --- a/src/SIM/SIMoptions.C +++ b/src/SIM/SIMoptions.C @@ -66,6 +66,8 @@ void SIMoptions::setLinearSolver (const std::string& eqsolver) solver = SystemMatrix::SPR; else if (eqsolver == "superlu") solver = SystemMatrix::SPARSE; + else if (eqsolver == "umfpack") + solver = SystemMatrix::UMFPACK; else if (eqsolver == "samg") solver = SystemMatrix::SAMG; else if (eqsolver == "petsc") @@ -259,6 +261,8 @@ bool SIMoptions::parseOldOptions (int argc, char** argv, int& i) } else if (!strcmp(argv[i],"-samg")) solver = SystemMatrix::SAMG; + else if (!strcmp(argv[i],"-umfpack")) + solver = SystemMatrix::UMFPACK; else if (!strcmp(argv[i],"-petsc")) solver = SystemMatrix::PETSC; else if (!strcmp(argv[i],"-istl"))