Tweaks to SuperLU constructors (#4063)

* Construct unwrapped Python object at Python construction time.

* std::move

* clang-format

* Don't use moved shared_ptr after move.

* Fix

* Update logging.

* Remove assert, reference.

* clang-format

* Tweak

* Do we really require pkgconfig at runtime? pkgconfig.py looks unused.

* Tweaks

Author: jhale

cpp/dolfinx/la/superlu_dist.cpp

@@ -70,7 +70,7 @@ template <typename T>
 std::unique_ptr<SuperLUDistStructs::SuperMatrix, SuperMatrixDeleter>
 create_supermatrix(const auto& A, auto& rowptr, auto& cols)
 {
-  spdlog::info("Start set_operator");
+  spdlog::info("Start create_supermatrix");
 
   auto map0 = A.index_map(0);
   auto map1 = A.index_map(1);
@@ -113,7 +113,7 @@ create_supermatrix(const auto& A, auto& rowptr, auto& cols)
   else
     static_assert(dependent_false_v<T>, "Invalid scalar type");
 
-  spdlog::info("Finished set_operator");
+  spdlog::info("Finished create_supermatrix");
   return p;
 }
 } // namespace
@@ -122,10 +122,12 @@ create_supermatrix(const auto& A, auto& rowptr, auto& cols)
 template <typename T>
 SuperLUDistMatrix<T>::SuperLUDistMatrix(std::shared_ptr<const MatrixCSR<T>> A,
                                         bool verbose)
-    : _matA(A),
-      _cols(std::make_unique<SuperLUDistStructs::vec_int_t>(col_indices(*A))),
-      _rowptr(std::make_unique<SuperLUDistStructs::vec_int_t>(row_indices(*A))),
-      _supermatrix(create_supermatrix<T>(*A, *_rowptr, *_cols)),
+    : _matA(std::move(A)),
+      _cols(
+          std::make_unique<SuperLUDistStructs::vec_int_t>(col_indices(*_matA))),
+      _rowptr(
+          std::make_unique<SuperLUDistStructs::vec_int_t>(row_indices(*_matA))),
+      _supermatrix(create_supermatrix<T>(*_matA, *_rowptr, *_cols)),
       _verbose(verbose)
 {
 }
@@ -137,6 +139,19 @@ SuperLUDistStructs::SuperMatrix* SuperLUDistMatrix<T>::supermatrix() const
   return _supermatrix.get();
 }
 
+//----------------------------------------------------------------------------
+template <typename T>
+const MatrixCSR<T>& SuperLUDistMatrix<T>::matA() const
+{
+  return *_matA;
+}
+
+//----------------------------------------------------------------------------
+template class la::SuperLUDistMatrix<double>;
+template class la::SuperLUDistMatrix<float>;
+template class la::SuperLUDistMatrix<std::complex<double>>;
+//----------------------------------------------------------------------------
+
 //----------------------------------------------------------------------------
 // Trick for declaring anonymous typedef structs from SuperLU_DIST
 struct dolfinx::la::SuperLUDistStructs::gridinfo_t : public ::gridinfo_t
@@ -153,11 +168,11 @@ void GridInfoDeleter::operator()(
 
 //----------------------------------------------------------------------------
 template <typename T>
-SuperLUDistSolver<T>::SuperLUDistSolver(std::shared_ptr<const MatrixCSR<T>> A,
-                                        bool verbose)
-    : _superlu_matA(SuperLUDistMatrix<T>(A, verbose)),
+SuperLUDistSolver<T>::SuperLUDistSolver(
+    std::shared_ptr<const SuperLUDistMatrix<T>> A, bool verbose)
+    : _superlu_matA(std::move(A)),
       _gridinfo(
-          [comm = A->comm()]
+          [comm = _superlu_matA->matA().comm()]
           {
             int nprow = dolfinx::MPI::size(comm);
             int npcol = 1;
@@ -174,8 +189,8 @@ SuperLUDistSolver<T>::SuperLUDistSolver(std::shared_ptr<const MatrixCSR<T>> A,
 template <typename T>
 int SuperLUDistSolver<T>::solve(const la::Vector<T>& b, la::Vector<T>& u) const
 {
-  int_t m = _superlu_matA.supermatrix()->nrow;
-  int_t m_loc = ((NRformat_loc*)(_superlu_matA.supermatrix()->Store))->m_loc;
+  int_t m = _superlu_matA->supermatrix()->nrow;
+  int_t m_loc = ((NRformat_loc*)(_superlu_matA->supermatrix()->Store))->m_loc;
 
   // RHS
   int_t ldb = m_loc;
@@ -207,7 +222,7 @@ int SuperLUDistSolver<T>::solve(const la::Vector<T>& b, la::Vector<T>& u) const
     dSOLVEstruct_t SOLVEstruct;
 
     spdlog::info("Call SuperLU_DIST pdgssvx()");
-    pdgssvx(&options, _superlu_matA.supermatrix(), &ScalePermstruct,
+    pdgssvx(&options, _superlu_matA->supermatrix(), &ScalePermstruct,
             u.array().data(), ldb, nrhs, _gridinfo.get(), &LUstruct,
             &SOLVEstruct, berr.data(), &stat, &info);
 
@@ -226,7 +241,7 @@ int SuperLUDistSolver<T>::solve(const la::Vector<T>& b, la::Vector<T>& u) const
     sSOLVEstruct_t SOLVEstruct;
 
     spdlog::info("Call SuperLU_DIST psgssvx()");
-    psgssvx(&options, _superlu_matA.supermatrix(), &ScalePermstruct,
+    psgssvx(&options, _superlu_matA->supermatrix(), &ScalePermstruct,
             u.array().data(), ldb, nrhs, _gridinfo.get(), &LUstruct,
             &SOLVEstruct, berr.data(), &stat, &info);
 
@@ -245,7 +260,7 @@ int SuperLUDistSolver<T>::solve(const la::Vector<T>& b, la::Vector<T>& u) const
     zSOLVEstruct_t SOLVEstruct;
 
     spdlog::info("Call SuperLU_DIST pzgssvx()");
-    pzgssvx(&options, _superlu_matA.supermatrix(), &ScalePermstruct,
+    pzgssvx(&options, _superlu_matA->supermatrix(), &ScalePermstruct,
             reinterpret_cast<doublecomplex*>(u.array().data()), ldb, nrhs,
             _gridinfo.get(), &LUstruct, &SOLVEstruct, berr.data(), &stat,
             &info);

cpp/dolfinx/la/superlu_dist.h

@@ -53,9 +53,12 @@ class SuperLUDistMatrix
   /// Copy assignment
   SuperLUDistMatrix& operator=(const SuperLUDistMatrix&) = delete;
 
-  /// Get non-const pointer to SuperLU_DIST native SuperMatrix.
+  /// Get non-const pointer to SuperLU_DIST SuperMatrix.
   SuperLUDistStructs::SuperMatrix* supermatrix() const;
 
+  /// Get MatrixCSR (const).
+  const MatrixCSR<T>& matA() const;
+
 private:
   // Saved matrix operator with rows and cols in required integer type.
   // cols and rowptr are required in opaque type "int_t" of
@@ -93,7 +96,7 @@ class SuperLUDistSolver
   /// @tparam T Scalar type.
   /// @param A Matrix to solve for.
   /// @param verbose Verbose output.
-  SuperLUDistSolver(std::shared_ptr<const MatrixCSR<T>> A,
+  SuperLUDistSolver(std::shared_ptr<const SuperLUDistMatrix<T>> A,
                     bool verbose = false);
 
   /// Copy constructor
@@ -113,7 +116,7 @@ class SuperLUDistSolver
 
 private:
   // Wrapped SuperLU SuperMatrix
-  const SuperLUDistMatrix<T> _superlu_matA;
+  std::shared_ptr<const SuperLUDistMatrix<T>> _superlu_matA;
 
   // Pointer to struct gridinfo_t
   std::unique_ptr<SuperLUDistStructs::gridinfo_t, GridInfoDeleter> _gridinfo;

python/doc/source/installation.rst

@@ -31,11 +31,11 @@ C++
 - `Boost <https://www.boost.org>`_
 - `CMake <https://cmake.org>`_ [build dependency]
 - HDF5 (with MPI support enabled)
-- MPI (MPI-3 or later).
-- `pkg-config <https://www.freedesktop.org/wiki/Software/pkg-config/>`_
+- MPI (MPI-3 or later)
+- `pkg-config <https://www.freedesktop.org/wiki/Software/pkg-config/>`_ [build dependency]
 - `pugixml <https://pugixml.org/>`_
 - `spdlog <https://github.com/gabime/spdlog/>`_
-- UFCx [``ufcx.h``, provided by FFCx package or FFCx UFCx CMake install
+- UFCx [``ufcx.h``, provided by FFCx Python package or FFCx UFCx CMake install
   at ``ffcx/cmake/*``]
 - At least one of ParMETIS [2]_, KaHIP or PT-SCOTCH [2]_
 
@@ -45,8 +45,9 @@ From ParMETIS, KaHIP or PT-SCOTCH, ParMETIS is recommended.
 
 - `ADIOS2 <https://github.com/ornladios/ADIOS2/>`_ (additional parallel
   IO support)
-- `PETSc <https://petsc.org/>`_ [1]_
-- `SLEPc <https://slepc.upv.es/>`_ (eigenvalue computations)
+- `PETSc <https://petsc.org/>`_ [1]_ (linear and non-linear problems)
+- `SLEPc <https://slepc.upv.es/>`_ (eigenvalue problems)
+- `SuperLU_DIST <https://github.com/xiaoyeli/superlu_dist/>`_ [2]_ (linear problems with ``dolfinx::la::MatrixCSR``).
 
 .. rubric:: Optional for demos
 
@@ -57,28 +58,43 @@ PETSc and FFCx are optional but recommended.
 Python interface
 ****************
 
-Requirements for the Python interface, in addition to the C++
-requirements.
+Requirements for the Python interface. Please see ``python/pyproject.toml`` for
+precise specification. Below we use the `pypi <https://pypi.org>`_ package names.
 
-.. rubric:: Required
+.. rubric:: Build system requirements
 
 - Python
-- Python CFFI (https://cffi.readthedocs.io/)
-- FFCx, UFL and Basix Python interface.
-- mpi4py (https://mpi4py.readthedocs.io/)
-- nanobind (https://github.com/wjakob/nanobind)
-- NumPy (https://www.numpy.org)
-- scikit-build-core[pyproject] (https://scikit-build-core.readthedocs.io)
+- DOLFINx C++ interface and all requirements
+- `scikit-build-core[pyproject] <https://scikit-build-core.readthedocs.io>`_
+- `mpi4py <https://mpi4py.readthedocs.io/>`_
+- `nanobind <https://github.com/wjakob/nanobind>`_ (static linking)
+- petsc4py (recommended, optional)
 
-.. rubric:: Optional
+.. rubric:: Required runtime dependencies
+
+- Python
+- fenics-basix, fenics-ffcx and fenics-ufl
+- `cffi <https://cffi.readthedocs.io/>`_
+- `mpi4py <https://mpi4py.readthedocs.io/>`_
+- `numpy <https://www.numpy.org>`_
+- `pkg-config <https://www.freedesktop.org/wiki/Software/pkg-config/>`_
+
+.. rubric:: Optional runtime dependencies
 
-- petsc4py (recommended)
+- petsc4py (linear and non-linear problems, recommended)
+- numba (custom kernels and assemblers)
+- `pyamg <https://github.com/pyamg/pyamg>`_ + scipy (serial linear problems)
 
 .. rubric:: Optional for demos
 
-- Numba
+- gmsh
+- networkx 
+- numba
+- matplotlib
+- petsc4py
 - pyamg
-- pyvista (for plotting)
+- pyvista
+- scipy
 - slepc4py
 
 Building and installing
@@ -118,4 +134,4 @@ the Python interface can be installed using::
        additionally configure MUMPS via PETSc with
        ``--download-mumps-avoid-mpi-in-place``.
 
-.. [2] PETSc can download and configure and build these libraries.
+.. [2] PETSc can also download, configure and build these libraries.

python/dolfinx/wrappers/la.cpp

@@ -229,7 +229,13 @@ void declare_superlu_dist_solver(nb::module_& m, const std::string& type)
           [](dolfinx::la::SuperLUDistSolver<T>* solver,
              std::shared_ptr<const dolfinx::la::MatrixCSR<T>> Amat,
              bool verbose)
-          { new (solver) dolfinx::la::SuperLUDistSolver<T>(Amat, verbose); },
+          {
+            auto A_superlu
+                = std::make_shared<const dolfinx::la::SuperLUDistMatrix<T>>(
+                    std::move(Amat), verbose);
+            new (solver) dolfinx::la::SuperLUDistSolver<T>(std::move(A_superlu),
+                                                           verbose);
+          },
           nb::arg("A"), nb::arg("verbose"))
       .def("solve", &dolfinx::la::SuperLUDistSolver<T>::solve);
 }