add umfpack solver to SparseMatrix

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)

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)
+if(SuiteSparse_UMFPACK_FOUND AND ISTL_FOUND)
-  find_package(SuiteSparse COMPONENTS umfpack)
-endif()
-
-if(SuiteSparse_UMFPACK_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)

src/IFEM.C

@@ -24,6 +24,9 @@
 #elif defined(HAVE_MPI)
 #include <mpi.h>
 #endif
+#ifdef HAS_UMFPACK
+#include <umfpack.h>
+#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

src/LinAlg/SparseMatrix.C

@@ -22,6 +22,9 @@
 #ifdef HAS_SAMG
 #include "samg.h"
 #endif
+#ifdef HAS_UMFPACK
+#include <umfpack.h>
+#endif
 #ifdef USE_OPENMP
 #include <omp.h>
 #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<IntVec>& 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++)

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;
 

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

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,

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"))