superlu_dist.h

// Copyright (C) 2026 Jack S. Hale, Chris N. Richardson
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#ifdef HAS_SUPERLU_DIST
 
#include <complex>
#include <dolfinx/common/MPI.h>
#include <dolfinx/la/MatrixCSR.h>
#include <dolfinx/la/Vector.h>
#include <memory>
#include <string>
#include <string_view>
#include <vector>
 
namespace dolfinx::la
{
class SuperLUDistStructs
{
public:
  struct SuperMatrix;
  struct vec_int_t;
  struct gridinfo_t;
  struct superlu_dist_options_t;
 
  struct sScalePermstruct_t;
  struct dScalePermstruct_t;
  struct zScalePermstruct_t;
 
  struct sLUstruct_t;
  struct dLUstruct_t;
  struct zLUstruct_t;
 
  struct sSOLVEstruct_t;
  struct dSOLVEstruct_t;
  struct zSOLVEstruct_t;
};
 
// SuperLU_DIST has structs that are 'typed' with prefixes d, s, z. This allows
// the solver class to select the typed set based on T.
namespace impl
{
template <typename...>
constexpr bool always_false_v = false;
 
template <typename T>
struct map
{
  static_assert(always_false_v<T>, "Invalid scalar type");
};

template <>
struct map<double>
{
  using ScalePermstruct_t = SuperLUDistStructs::dScalePermstruct_t;
  using LUstruct_t = SuperLUDistStructs::dLUstruct_t;
  using SOLVEstruct_t = SuperLUDistStructs::dSOLVEstruct_t;
};

template <>
struct map<float>
{
  using ScalePermstruct_t = SuperLUDistStructs::sScalePermstruct_t;
  using LUstruct_t = SuperLUDistStructs::sLUstruct_t;
  using SOLVEstruct_t = SuperLUDistStructs::sSOLVEstruct_t;
};

template <>
struct map<std::complex<double>>
{
  using ScalePermstruct_t = SuperLUDistStructs::zScalePermstruct_t;
  using LUstruct_t = SuperLUDistStructs::zLUstruct_t;
  using SOLVEstruct_t = SuperLUDistStructs::zSOLVEstruct_t;
};
 
} // namespace impl

template <typename T>
using map_t = impl::map<T>;

struct SuperMatrixDeleter
{
  void operator()(SuperLUDistStructs::SuperMatrix* A) const noexcept;
};

template <typename T>
class SuperLUDistMatrix
{
public:
  SuperLUDistMatrix(const MatrixCSR<T>& A);

  SuperLUDistMatrix(const SuperLUDistMatrix&) = delete;

  SuperLUDistMatrix& operator=(const SuperLUDistMatrix&) = delete;

  MPI_Comm comm() const;

  SuperLUDistStructs::SuperMatrix* supermatrix() const;
 
private:
  dolfinx::MPI::Comm _comm;
  // Deep copy of values from MatrixCSR.
  std::vector<T> _matA_values;
  // cols and rowptr are required in opaque type "int_t" of
  // SuperLU_DIST.
  std::unique_ptr<SuperLUDistStructs::vec_int_t> _cols;
  std::unique_ptr<SuperLUDistStructs::vec_int_t> _rowptr;
 
  // Pointer to native SuperMatrix for use in solver
  std::unique_ptr<SuperLUDistStructs::SuperMatrix, SuperMatrixDeleter>
      _supermatrix;
};

struct GridInfoDeleter
{
  void operator()(SuperLUDistStructs::gridinfo_t* g) const noexcept;
};

struct ScalePermStructDeleter
{
  void operator()(SuperLUDistStructs::dScalePermstruct_t* s) const noexcept;
  void operator()(SuperLUDistStructs::sScalePermstruct_t* s) const noexcept;
  void operator()(SuperLUDistStructs::zScalePermstruct_t* s) const noexcept;
};

struct LUStructDeleter
{
  void operator()(SuperLUDistStructs::dLUstruct_t* l) const noexcept;
  void operator()(SuperLUDistStructs::sLUstruct_t* l) const noexcept;
  void operator()(SuperLUDistStructs::zLUstruct_t* l) const noexcept;
};

struct SolveStructDeleter
{
  SuperLUDistStructs::superlu_dist_options_t* o;

  void operator()(SuperLUDistStructs::dSOLVEstruct_t* S) const noexcept;
  void operator()(SuperLUDistStructs::sSOLVEstruct_t* S) const noexcept;
  void operator()(SuperLUDistStructs::zSOLVEstruct_t* S) const noexcept;
};

template <typename T>
class SuperLUDistSolver
{
public:
  SuperLUDistSolver(std::shared_ptr<const SuperLUDistMatrix<T>> A);

  SuperLUDistSolver(const SuperLUDistSolver&) = delete;

  SuperLUDistSolver& operator=(const SuperLUDistSolver&) = delete;

  ~SuperLUDistSolver();

  void set_option(std::string_view name, std::string_view value);

  void set_options(const SuperLUDistStructs::superlu_dist_options_t& options);

  void set_A(std::shared_ptr<const SuperLUDistMatrix<T>> A,
             std::string_view fact);

  int solve(const Vector<T>& b, Vector<T>& u);
 
private:
  // Assembled left-hand side matrix
  std::shared_ptr<const SuperLUDistMatrix<T>> _superlu_matA;
 
  // Pointer to struct superlu_dist_options_t
  std::unique_ptr<SuperLUDistStructs::superlu_dist_options_t> _options;
 
  // Pointer to struct gridinfo_t
  std::unique_ptr<SuperLUDistStructs::gridinfo_t, GridInfoDeleter> _gridinfo;
 
  // Pointer to 'typed' struct *ScalePermstruct_t
  std::unique_ptr<typename map_t<T>::ScalePermstruct_t, ScalePermStructDeleter>
      _scalepermstruct;
  // Pointer to 'typed' struct *LUstruct_t
  std::unique_ptr<typename map_t<T>::LUstruct_t, LUStructDeleter> _lustruct;
  // Pointer to 'typed' struct *SOLVEstruct
  std::unique_ptr<typename map_t<T>::SOLVEstruct_t, SolveStructDeleter>
      _solvestruct;
 
  // True once pdgssvx has populated LUstruct with per-block-column arrays
  // that must be released via Destroy_LU before LUstructFree.
  bool _factored = false;
};
} // namespace dolfinx::la
#endif