petsc.h

// Copyright (C) 2004-2026 Johan Hoffman, Johan Jansson, Anders Logg,
// Garth N. Wells and Jack S. Hale
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#ifdef HAS_PETSC
 
#include "Vector.h"
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <dolfinx/common/petsc.h>
#include <dolfinx/common/types.h>
#include <functional>
#include <iterator>
#include <optional>
#include <petscksp.h>
#include <petscmat.h>
#include <petscvec.h>
#include <ranges>
#include <span>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
 
namespace dolfinx::common
{
class IndexMap;
} // namespace dolfinx::common
 
namespace dolfinx::la
{
class SparsityPattern;

namespace petsc
{
std::vector<Vec>
create_vectors(MPI_Comm comm,
               const std::vector<std::span<const PetscScalar>>& x);

Vec create_vector(const common::IndexMap& map, int bs);

Vec create_vector(MPI_Comm comm, std::array<std::int64_t, 2> range,
                  std::span<const std::int64_t> ghosts, int bs);

Vec create_vector_wrap(const common::IndexMap& map, int bs,
                       std::span<const PetscScalar> x);

template <class V>
Vec create_vector_wrap(const la::Vector<V>& x)
{
  assert(x.index_map());
  return create_vector_wrap(*x.index_map(), x.bs(), x.array());
}

std::vector<IS> create_index_sets(
    const std::vector<
        std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);

std::vector<IS> create_global_index_sets(
    const std::vector<
        std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);

std::vector<std::vector<PetscScalar>> get_local_vectors(
    const Vec x,
    const std::vector<
        std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);

void scatter_local_vectors(
    Vec x, const std::vector<std::span<const PetscScalar>>& x_b,
    const std::vector<
        std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);

Mat create_matrix(MPI_Comm comm, const SparsityPattern& sp,
                  std::optional<std::string_view> type = std::nullopt,
                  std::optional<ISLocalToGlobalMapping> rlgmap = std::nullopt,
                  std::optional<ISLocalToGlobalMapping> clgmap = std::nullopt);

MatNullSpace create_nullspace(MPI_Comm comm, std::span<const Vec> basis);

class Vector
{
public:
  Vector(const common::IndexMap& map, int bs);
 
  // Delete copy constructor to avoid accidental copying of 'heavy' data
  Vector(const Vector& x) = delete;

  Vector(Vector&& x) noexcept;

  Vector(Vec x, bool inc_ref_count);

  ~Vector();
 
  // Assignment operator (disabled)
  Vector& operator=(const Vector& x) = delete;

  Vector& operator=(Vector&& x) noexcept;

  Vector copy() const;

  std::int64_t size() const;

  std::int32_t local_size() const;

  std::array<std::int64_t, 2> local_range() const;

  MPI_Comm comm() const;

  void set_options_prefix(std::string_view options_prefix);

  std::string get_options_prefix() const;

  void set_from_options();

  Vec vec() const;
 
private:
  // PETSc Vec pointer
  Vec _x;
};

class Matrix
{
public:
  static auto set_fn(Mat A, InsertMode mode)
  {
    return [A, mode, cache = std::vector<PetscInt>()](
               const common::LocalIndexRange auto& rows,
               const common::LocalIndexRange auto& cols,
               std::span<const PetscScalar> vals) mutable -> int
    {
      PetscErrorCode ierr;
#ifdef PETSC_USE_64BIT_INDICES
      cache.resize(std::ranges::size(rows) + std::ranges::size(cols));
      std::ranges::copy(rows, cache.begin());
      std::ranges::copy(cols,
                        std::next(cache.begin(), std::ranges::size(rows)));
      const PetscInt* _rows = cache.data();
      const PetscInt* _cols = cache.data() + std::ranges::size(rows);
      ierr = MatSetValuesLocal(A, std::ranges::size(rows), _rows,
                               std::ranges::size(cols), _cols, vals.data(),
                               mode);
#else
      ierr = MatSetValuesLocal(A, std::ranges::size(rows),
                               std::ranges::data(rows), std::ranges::size(cols),
                               std::ranges::data(cols), vals.data(), mode);
#endif
 
#ifndef NDEBUG
      common::petsc::check(ierr, "MatSetValuesLocal");
#endif
      return ierr;
    };
  }

  static auto set_block_fn(Mat A, InsertMode mode)
  {
    return [A, mode, cache = std::vector<PetscInt>()](
               const common::LocalIndexRange auto& rows,
               const common::LocalIndexRange auto& cols,
               std::span<const PetscScalar> vals) mutable -> int
    {
      PetscErrorCode ierr;
#ifdef PETSC_USE_64BIT_INDICES
      cache.resize(std::ranges::size(rows) + std::ranges::size(cols));
      std::ranges::copy(rows, cache.begin());
      std::ranges::copy(cols,
                        std::next(cache.begin(), std::ranges::size(rows)));
      const PetscInt* _rows = cache.data();
      const PetscInt* _cols = cache.data() + std::ranges::size(rows);
      ierr = MatSetValuesBlockedLocal(A, std::ranges::size(rows), _rows,
                                      std::ranges::size(cols), _cols,
                                      vals.data(), mode);
#else
      ierr = MatSetValuesBlockedLocal(
          A, std::ranges::size(rows), std::ranges::data(rows),
          std::ranges::size(cols), std::ranges::data(cols), vals.data(), mode);
#endif
 
#ifndef NDEBUG
      common::petsc::check(ierr, "MatSetValuesBlockedLocal");
#endif
      return ierr;
    };
  }

  static auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)
  {
    return [A, bs0, bs1, mode, cache0 = std::vector<PetscInt>(),
            cache1 = std::vector<PetscInt>()](
               const common::LocalIndexRange auto& rows,
               const common::LocalIndexRange auto& cols,
               std::span<const PetscScalar> vals) mutable -> int
    {
      PetscErrorCode ierr;
      const std::int32_t* _rows = std::ranges::data(rows);
      const std::int32_t* _cols = std::ranges::data(cols);
      const std::size_t num_rows = std::ranges::size(rows);
      const std::size_t num_cols = std::ranges::size(cols);
      cache0.resize(bs0 * num_rows);
      cache1.resize(bs1 * num_cols);
      for (std::size_t i = 0; i < num_rows; ++i)
        for (int k = 0; k < bs0; ++k)
          cache0[bs0 * i + k] = bs0 * _rows[i] + k;
 
      for (std::size_t i = 0; i < num_cols; ++i)
        for (int k = 0; k < bs1; ++k)
          cache1[bs1 * i + k] = bs1 * _cols[i] + k;
 
      ierr = MatSetValuesLocal(A, cache0.size(), cache0.data(), cache1.size(),
                               cache1.data(), vals.data(), mode);
#ifndef NDEBUG
      common::petsc::check(ierr, "MatSetValuesLocal");
#endif
      return ierr;
    };
  }

  Matrix(MPI_Comm comm, const SparsityPattern& sp,
         std::optional<std::string_view> type = std::nullopt);

  Matrix(Mat A, bool inc_ref_count);
 
  // Copy constructor (deleted)
  Matrix(const Matrix& A) = delete;

  Matrix(Matrix&& A) noexcept;

  ~Matrix();
 
  // Assignment operator (deleted)
  Matrix& operator=(const Matrix& A) = delete;

  Matrix& operator=(Matrix&& A) noexcept;

  std::array<std::int64_t, 2> size() const;

  Vec create_vector(std::size_t dim) const;

  Mat mat() const;
 
  //--- Special PETSc Functions ---

  void set_options_prefix(std::string_view options_prefix);

  std::string get_options_prefix() const;

  void set_from_options();
 
private:
  // PETSc Mat pointer
  Mat _matA;
};

class KrylovSolver
{
public:
  explicit KrylovSolver(MPI_Comm comm);

  KrylovSolver(KSP ksp, bool inc_ref_count);
 
  // Copy constructor (deleted)
  KrylovSolver(const KrylovSolver& solver) = delete;

  KrylovSolver(KrylovSolver&& solver) noexcept;

  ~KrylovSolver();
 
  // Assignment operator (deleted)
  KrylovSolver& operator=(const KrylovSolver&) = delete;

  KrylovSolver& operator=(KrylovSolver&& solver) noexcept;

  void set_operator(const Mat A);

  void set_operators(const Mat A, const Mat P);

  [[nodiscard(
      "check the converged reason - positive on convergence, negative on "
      "divergence")]] KSPConvergedReason
  solve(Vec x, const Vec b, bool transpose = false);

  void set_options_prefix(std::string_view options_prefix);

  std::string get_options_prefix() const;

  void set_from_options() const;

  KSP ksp() const;
 
private:
  // PETSc solver pointer
  KSP _ksp;
};
} // namespace petsc
} // namespace dolfinx::la
 
#endif