Linear algebra (dolfinx::la)#

namespace la#

Linear algebra interface.

Interface to linear algebra data structures and solvers.

Typedefs

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

Map scalar type to SuperLU_DIST ‘typed’ structs.

Enums

enum class BlockMode : int#

Modes for representing block structured matrices.

Values:

enumerator compact#
enumerator expanded#

Each entry in the sparsity pattern of the matrix refers to a block of data of size (bs[0], bs[1]).

enum class Norm : std::int8_t#

Norm types.

Values:

enumerator l1#
enumerator l2#
enumerator linf#
enumerator frobenius#

Functions

template<typename T>
dolfinx::la::MatrixCSR<T> matmul(const dolfinx::la::MatrixCSR<T> &A, const dolfinx::la::MatrixCSR<T> &B)#

Compute C = A*B as a distributed MatrixCSR.

Note

Currently only supports block-size 1 matrices.

Parameters:

  • A – Left matrix.

  • B – Right matrix.

Returns:

The product C = A*B as a MatrixCSR with row distribution matching A and column distribution determined by B. The row IndexMap of C has no ghosts, and the column IndexMap of C will generally be a larger superset of the column IndexMap of B.

template<typename T, int BS0 = -1, int BS1 = -1>
dolfinx::la::MatrixCSR<T> transpose(const dolfinx::la::MatrixCSR<T> &A)#

Compute the distributed transpose of a MatrixCSR.

Given A with row map R (n_row owned rows) and column map C (n_col owned columns plus any ghost columns from remote ranks), returns Aᵀ with: row map = ghost-free column map of A (n_col owned rows, no ghosts) col map = row map of A extended with any ghost row indices that appear as nonzeros contributed by remote ranks.

Communication strategy:

  1. Exchange per-neighbour entry counts (MPI_Neighbor_alltoall).

  2. Build displacements and pack send buffers.

  3. Post non-blocking data exchange (three MPI_Ineighbor_alltoallv in flight simultaneously: row global indices, column global indices, and values).

  4. While data exchange is in flight, compute the local (diagonal-block) part of the transpose via local_transpose.

  5. Wait for data, merge received entries, finalise sparsity and construct the MatrixCSR.

Parameters:

A – Input matrix.

Template Parameters:

  • T – Scalar type

  • BS0 – Row block size of A, or -1 for unoptimized.

  • BS1 – Col block size of A, or -1 for unoptimized.

Returns:

Aᵀ as a distributed MatrixCSR.

template<class V>
auto inner_product(const V &a, const V &b)#

Compute the inner product of two vectors.

Computes a^{H} b (a^{T} b if a and b are real). The two vectors must have the same parallel layout.

Note

Collective MPI operation

Parameters:
Returns:

Inner product between a and b.

template<class V>
auto squared_norm(const V &a)#

Compute the squared L2 norm of vector.

Note

Collective MPI operation.

template<class V>
auto norm(const V &x, Norm type = Norm::l2)#

Compute the norm of the vector.

Note

Collective MPI operation.

Parameters:

  • x – Vector to compute the norm of.

  • type – Norm type.

template<class V>
void orthonormalize(std::vector<std::reference_wrapper<V>> basis)#

Orthonormalize a set of vectors.

Template Parameters:

V – dolfinx::la::Vector

Parameters:

basis – [inout] The set of vectors to orthonormalise. The vectors must have identical parallel layouts. The vectors are modified in-place.

template<class V>
bool is_orthonormal(std::vector<std::reference_wrapper<const V>> basis, dolfinx::scalar_value_t<typename V::value_type> eps = std::numeric_limits<dolfinx::scalar_value_t<typename V::value_type>>::epsilon())#

Test if basis is orthonormal.

Returns true if ||x_i - x_j|| - delta_{ij} < eps for all i, j, and otherwise false.

Parameters:

  • basis – [in] Set of vectors to check.

  • eps – [in] Tolerance.

Returns:

True is basis is orthonormal, otherwise false.

struct GridInfoDeleter#
#include <superlu_dist.h>

Call library cleanup and delete pointer. For use with std::unique_ptr holding gridinfo_t.

Public Functions

void operator()(SuperLUDistStructs::gridinfo_t *g) const noexcept#

Deletion of gridinfo_t.

Parameters:

g –

struct LUStructDeleter#
#include <superlu_dist.h>

Call library cleanup and delete pointer. For use with std::unique_ptr holding *LUstruct_t

Public Functions

void operator()(SuperLUDistStructs::dLUstruct_t *l) const noexcept#

double implementation

void operator()(SuperLUDistStructs::sLUstruct_t *l) const noexcept#

float implementation

void operator()(SuperLUDistStructs::zLUstruct_t *l) const noexcept#

doublecomplex implementation

template<typename Scalar, typename Container = std::vector<Scalar>, typename ColContainer = std::vector<std::int32_t>, typename RowPtrContainer = std::vector<std::int64_t>>
class MatrixCSR#
#include <MatrixCSR.h>

Distributed sparse matrix using compressed sparse row storage.

The matrix storage format is compressed sparse row, and is partitioned by row across MPI ranks. The matrix can be assembled into using the usual DOLFINx assembly routines. Matrix internal data can be accessed for interfacing with other code.

Warning

The class is experimental and subject to change.

Template Parameters:

  • Scalar – Scalar type of matrix entries

  • Container – Container type for storing matrix entries

  • ColContainer – Column index container type

  • RowPtrContainer – Row pointer container type

Public Types

using value_type = Scalar#

Scalar type.

using container_type = Container#

Matrix entries container type.

using column_container_type = ColContainer#

Column index container type.

using rowptr_container_type = RowPtrContainer#

Row pointer container type.

Public Functions

template<int BS0 = 1, int BS1 = 1>
inline auto mat_set_values()#

Insertion functor for setting values in a matrix. It is typically used in finite element assembly functions.

Create a function to set values in a MatrixCSR. The function signature is int mat_set_fn(std::span<const std::int32_t rows, std::span<const std::int32_t cols, std::span<const value_type> data). The rows and columns use process local indexing, and the given rows and columns must pre-exist in the sparsity pattern of the matrix. Insertion into “ghost” rows (in the ghost region of the row IndexMap) is permitted, so long as there are correct entries in the sparsity pattern.

Note

Using rows or columns which are not in the sparsity will result in undefined behaviour (or an assert failure in Debug mode).

Note

Matrix block size may be (1, 1) or (BS0, BS1)

Note

Data block size may be (1, 1) or (BS0, BS1)

Template Parameters:

  • BS0 – Row block size of data for insertion

  • BS1 – Column block size of data for insertion

Returns:

Function for inserting values into A

template<int BS0 = 1, int BS1 = 1>
inline auto mat_add_values()#

Insertion functor for adding values to a matrix. It is typically used in finite element assembly functions.

Create a function to add values to a MatrixCSR. The function signature is int mat_add_fn(std::span<const std::int32_t rows, std::span<const std::int32_t cols, std::span<const value_type> data). The rows and columns use process local indexing, and the given rows and columns must pre-exist in the sparsity pattern of the matrix. Insertion into “ghost” rows (in the ghost region of the row IndexMap) is permitted, so long as there are correct entries in the sparsity pattern.

Note

Using rows or columns which are not in the sparsity will result in undefined behaviour (or an assert failure in Debug mode).

Note

Matrix block size may be (1, 1) or (BS0, BS1)

Note

Data block size may be (1, 1) or (BS0, BS1)

Template Parameters:

  • BS0 – Row block size of data for insertion

  • BS1 – Column block size of data for insertion

Returns:

Function for inserting values into A

template<SparsityImplementation T>
MatrixCSR(const T &p, BlockMode mode = BlockMode::compact)#

Create a distributed matrix.

The structure of the matrix depends entirely on the input SparsityPattern, which must be finalized. The matrix storage is distributed Compressed Sparse Row: the matrix is distributed by row across processes, and on each process, there is a list of column indices and matrix entries for each row stored. This exactly matches the layout of the SparsityPattern. There is some overlap of matrix rows between processes to allow for independent Finite Element assembly, after which, the ghost rows should be sent to the row owning processes by calling scatter_rev().

Note

The block size of the matrix is given by the block size of the input SparsityPattern.

Parameters:

  • p – [in] The sparsity pattern which describes the parallel distribution and the non-zero structure.

  • mode – [in] Block mode. When the block size is greater than one, the storage can be “compact” where each matrix entry refers to a block of data (stored row major), or “expanded” where each matrix entry is individual. In the “expanded” case, the sparsity is expanded for every entry in the block, and the block size of the matrix is set to (1, 1).

MatrixCSR(MatrixCSR &&A) = default#

Move constructor

Todo:

Check handling of MPI_Request

MatrixCSR(const MatrixCSR &A) = default#

Copy constructor

Todo:

Check handling of MPI_Request

template<typename Scalar0, typename Container0, typename ColContainer0, typename RowPtrContainer0>
inline explicit MatrixCSR(const MatrixCSR<Scalar0, Container0, ColContainer0, RowPtrContainer0> &A)#

Copy-convert matrix, possibly using to different container types.

Examples of use for this constructor include copying a matrix to a different value type, or copying the matrix to a GPU.

Template Parameters:

Mat –

Parameters:

A – Matrix to copy.

inline void set(value_type x)#

Set all non-zero local entries to a value, including entries in ghost rows.

Deprecated:

Use std::ranges::fill(A.values(), x) instead.

Parameters:

x – [in] The value to set non-zero matrix entries to

template<int BS0, int BS1>
inline void set(std::span<const value_type> x, std::span<const std::int32_t> rows, std::span<const std::int32_t> cols)#

Set values in the matrix.

Note

Only entries included in the sparsity pattern used to initialize the matrix can be set.

Note

All indices are local to the calling MPI rank and entries cannot be set in ghost rows.

Note

This should be called after scatter_rev. Using before scatter_rev will set the values correctly, but incoming values may get added to them during a subsequent reverse scatter operation.

Template Parameters:

  • BS0 – Data row block size

  • BS1 – Data column block size

Parameters:

  • x – [in] The m by n dense block of values (row-major) to set in the matrix

  • rows – [in] The row indices of x

  • cols – [in] The column indices of x

template<int BS0 = 1, int BS1 = 1>
inline void add(std::span<const value_type> x, std::span<const std::int32_t> rows, std::span<const std::int32_t> cols)#

Accumulate values in the matrix.

Note

Only entries included in the sparsity pattern used to initialize the matrix can be accumulated into.

Note

All indices are local to the calling MPI rank and entries may go into ghost rows.

Note

Use scatter_rev after all entries have been added to send ghost rows to owners. Adding more entries after scatter_rev is allowed, but another call to scatter_rev will then be required.

Template Parameters:

  • BS0 – Row block size of data

  • BS1 – Column block size of data

Parameters:

  • x – [in] The m by n dense block of values (row-major) to add to the matrix

  • rows – [in] The row indices of x

  • cols – [in] The column indices of x

inline std::int32_t num_owned_rows() const#

Number of local rows excluding ghost rows.

inline std::int32_t num_all_rows() const#

Number of local rows including ghost rows.

inline std::vector<value_type> to_dense() const#

Copy to a dense matrix.

Note

This function is typically used for debugging and not used in production.

Note

Ghost rows are also returned, and these can be truncated manually by using num_owned_rows() if required.

Note

If the block size is greater than 1, the entries are expanded.

Returns:

Dense copy of the part of the matrix on the calling rank. Storage is row-major.

inline void scatter_rev()#

Transfer ghost row data to the owning ranks accumulating received values on the owned rows, and zeroing any existing data in ghost rows.

This process is analogous to scatter_rev for Vector except that the values are always accumulated on the owning process.

inline void scatter_rev_begin()#

Begin transfer of ghost row data to owning ranks, where it will be accumulated into existing owned rows.

Note

Calls to this function must be followed by MatrixCSR::scatter_rev_end(). Between the two calls matrix values must not be changed.

Note

This function does not change the matrix data. Data update only occurs with scatter_rev_end().

inline void scatter_rev_end()#

End transfer of ghost row data to owning ranks.

Note

Must be preceded by MatrixCSR::scatter_rev_begin().

Note

Matrix data received from other processes will be accumulated into locally owned rows, and ghost rows will be zeroed.

inline double squared_norm() const#

Compute the Frobenius norm squared across all processes.

Note

MPI Collective

void mult(Vector<value_type> &x, Vector<value_type> &y) const#

Compute the product y += Ax.

The vectors x and y must have parallel layouts that are compatible with A. In detail, x must have the same IndexMap as the matrix columns, A.index_map(1) and y must have the same owned indices as the matrix rows in A.index_map(0). Only owned entries of y are updated, so any ghost entries of y are not affected.

Computes y += A*x for a parallel CSR matrix A and parallel dense vectors x,y

Parameters:

  • x – [in] Vector to apply A to.

  • y – [inout] Vector to accumulate the result into.

void multT(Vector<value_type> &x, Vector<value_type> &y) const#

Compute the product y += A^T x.

Out-of-line implementation of MatrixCSR::multT.

Performs the distributed sparse matrix–vector product with the transpose of the matrix. The computation is split into two phases:

  1. Off-diagonal phase — contributions from the off-diagonal block (ghost columns of A, i.e. columns owned by remote ranks) are accumulated into the ghost region of y. The ghost entries of y are zeroed before this phase so that stale values do not pollute the result. A reverse scatter (scatter_rev) then reduces those ghost contributions back to the owning ranks.

  2. Diagonal phase — contributions from the diagonal block (locally owned columns of A) are accumulated into the owned entries of y.

Layout requirements:

  • x must share the same IndexMap as the matrix rows, A.index_map(0).

  • y must share the same owned indices as the matrix columns, A.index_map(1). Only owned entries of y are meaningful after the call; ghost entries are used as scratch and left in an unspecified state.

See also

MatrixCSR::multT for the full specification.

Note

y is accumulated into (not overwritten) on the owned entries. Zero y before the first call if a fresh result is required.

Parameters:

  • x – [in] Vector to apply A^T to; must cover the row space of A.

  • y – [inout] Vector accumulated into; covers the column space of A.

inline MPI_Comm comm() const#

Get MPI communicator that matrix is defined on.

inline std::shared_ptr<const common::IndexMap> index_map(int dim) const#

Index map for the row or column space.

The row index map contains ghost entries for rows which may be inserted into and the column index map contains all local and ghost columns that may exist in the owned rows.

Returns:

Row (0) or column (1) index map.

inline container_type &values()#

Get local data values.

Note

Includes ghost values.

inline const container_type &values() const#

Get local values (const version).

Note

Includes ghost values

inline const rowptr_container_type &row_ptr() const#

Get local row pointers.

Note

Includes pointers to ghost rows

inline const column_container_type &cols() const#

Get local column indices

Note

Includes columns in ghost rows

inline const rowptr_container_type &off_diag_offset() const#

Get the start of off-diagonal (unowned columns) on each row, allowing the matrix to be split (virtually) into two parts.

Operations (such as matrix-vector multiply) between the owned parts of the matrix and vector can then be performed separately from operations on the unowned parts.

Note

Includes ghost rows, which should be ‘truncated’ by the caller if not required.

inline std::array<int, 2> block_size() const#

Get block sizes.

Returns:

Block sizes for rows (0) and columns (1).

inline BlockMode block_mode() const#

Get ‘block mode’.

template<SparsityImplementation SparsityType>
MatrixCSR(const SparsityType &p, BlockMode mode)#

Create a distributed matrix.

The structure of the matrix depends entirely on the input SparsityPattern, which must be finalized. The matrix storage is distributed Compressed Sparse Row: the matrix is distributed by row across processes, and on each process, there is a list of column indices and matrix entries for each row stored. This exactly matches the layout of the SparsityPattern. There is some overlap of matrix rows between processes to allow for independent Finite Element assembly, after which, the ghost rows should be sent to the row owning processes by calling scatter_rev().

Note

The block size of the matrix is given by the block size of the input SparsityPattern.

Parameters:

  • p – [in] The sparsity pattern which describes the parallel distribution and the non-zero structure.

  • mode – [in] Block mode. When the block size is greater than one, the storage can be “compact” where each matrix entry refers to a block of data (stored row major), or “expanded” where each matrix entry is individual. In the “expanded” case, the sparsity is expanded for every entry in the block, and the block size of the matrix is set to (1, 1).

struct ScalePermStructDeleter#
#include <superlu_dist.h>

Call library cleanup and delete pointer. For use with std::unique_ptr holding *ScalePermstruct_t

Public Functions

void operator()(SuperLUDistStructs::dScalePermstruct_t *s) const noexcept#

double implementation

void operator()(SuperLUDistStructs::sScalePermstruct_t *s) const noexcept#

float implementation

void operator()(SuperLUDistStructs::zScalePermstruct_t *s) const noexcept#

doublecomplex implementation

class SLEPcEigenSolver#
#include <slepc.h>

This class provides an eigenvalue solver for PETSc matrices. It is a wrapper for the SLEPc eigenvalue solver.

Public Functions

explicit SLEPcEigenSolver(MPI_Comm comm)#

Create eigenvalue solver.

SLEPcEigenSolver(EPS eps, bool inc_ref_count)#

Create eigenvalue solver from EPS object.

SLEPcEigenSolver(SLEPcEigenSolver &&solver) noexcept#

Move constructor.

~SLEPcEigenSolver()#

Destructor.

SLEPcEigenSolver &operator=(SLEPcEigenSolver &&solver) noexcept#

Move assignment.

void set_operators(const Mat A, const Mat B)#

Set operators (B may be nullptr for regular eigenvalues problems)

void solve()#

Compute all eigenpairs of the matrix A (solve \(A x = \lambda x\)).

void solve(std::int64_t n)#

Compute the n first eigenpairs of the matrix A (solve \(A x = \lambda x\))

std::complex<PetscReal> get_eigenvalue(int i) const#

Get ith eigenvalue.

void get_eigenpair(PetscScalar &lr, PetscScalar &lc, Vec r, Vec c, int i) const#

Get ith eigenpair.

int get_iteration_number() const#

Get the number of iterations used by the solver.

std::int64_t get_number_converged() const#

Get the number of converged eigenvalues.

void set_options_prefix(const std::string &options_prefix)#

Sets the prefix used by PETSc when searching the PETSc options database

std::string get_options_prefix() const#

Returns the prefix used by PETSc when searching the PETSc options database

void set_from_options() const#

Set options from PETSc options database.

EPS eps() const#

Return SLEPc EPS pointer.

MPI_Comm comm() const#

Return MPI communicator.

struct SolveStructDeleter#
#include <superlu_dist.h>

Call library cleanup and delete pointer. For use with std::unique_ptr holding *SOLVEstruct_t

Public Functions

void operator()(SuperLUDistStructs::dSOLVEstruct_t *S) const noexcept#

double implementation

void operator()(SuperLUDistStructs::sSOLVEstruct_t *S) const noexcept#

float implementation

void operator()(SuperLUDistStructs::zSOLVEstruct_t *S) const noexcept#

doublecomplex implementation

Public Members

SuperLUDistStructs::superlu_dist_options_t *o#

Pointer to options - required for *SOLVEstruct_t cleanup function.

class SparsityPattern#
#include <SparsityPattern.h>

Sparsity pattern data structure that can be used to initialize sparse matrices. After assembly, column indices are always sorted in increasing order. Ghost entries are kept after assembly.

Public Functions

SparsityPattern(MPI_Comm comm, std::array<std::shared_ptr<const common::IndexMap>, 2> maps, std::array<int, 2> bs)#

Create an empty sparsity pattern with specified dimensions.

Parameters:

  • comm – [in] Communicator that the pattern is defined on.

  • maps – [in] Index maps describing the [0] row and [1] column index ranges (up to a block size).

  • bs – [in] Block sizes for the [0] row and [1] column maps.

SparsityPattern(MPI_Comm comm, const std::vector<std::vector<const SparsityPattern*>> &patterns, const std::array<std::vector<std::pair<std::reference_wrapper<const common::IndexMap>, int>>, 2> &maps, const std::array<std::vector<int>, 2> &bs)#

Create a new sparsity pattern by concatenating sub-patterns, e.g. pattern =[ pattern00 ][ pattern 01] [ pattern10 ][ pattern 11]

Parameters:

  • comm – [in] Communicator that the pattern is defined on.

  • patterns – [in] Rectangular array of sparsity pattern. The patterns must not be finalised. Null block are permitted/

  • maps – [in] Pairs of (index map, block size) for each row block (maps[0]) and column blocks (maps[1])/

  • bs – [in] Block sizes for the sparsity pattern entries/

SparsityPattern(SparsityPattern &&pattern) = default#

Move constructor.

~SparsityPattern() = default#

Destructor.

SparsityPattern &operator=(SparsityPattern &&pattern) = default#

Move assignment.

void insert(std::int32_t row, std::int32_t col)#

Insert non-zero locations using local (process-wise) indices.

Parameters:

  • row – [in] local row index

  • col – [in] local column index

void insert(std::span<const std::int32_t> rows, std::span<const std::int32_t> cols)#

Insert non-zero locations using local (process-wise) indices.

This routine inserts non-zero locations at the outer product of rows and cols into the sparsity pattern, i.e. adds the matrix entries at A[row[i], col[j]] for all i, j.

Parameters:

  • rows – [in] list of the local row indices

  • cols – [in] list of the local column indices

void insert_diagonal(std::span<const std::int32_t> rows)#

Insert non-zero locations on the diagonal.

Parameters:

rows – [in] Rows in local (process-wise) indices. The indices must exist in the row IndexMap.

void finalize()#

Finalize sparsity pattern and communicate off-process entries.

std::shared_ptr<const common::IndexMap> index_map(int dim) const#

Index map for given dimension dimension. Returns the index map for rows and columns that will be set by the current MPI rank.

Note

After finalization, the column index map is updated to account for additional column entries from other processes.

Parameters:

dim – [in] Requested map, row (0) or column (1).

Returns:

The index map.

std::vector<std::int64_t> column_indices() const#

Global indices of non-zero columns on owned rows.

Note

The ghosts are computed only once SparsityPattern::finalize has been called.

Returns:

Global index non-zero columns on this process, including ghosts.

int block_size(int dim) const#

Return index map block size for dimension dim.

std::int64_t num_nonzeros() const#

Number of nonzeros on this rank after assembly, including ghost rows.

std::int32_t nnz_diag(std::int32_t row) const#

Number of non-zeros in owned columns (diagonal block) on a given row.

Note

Can also be used on ghost rows

std::int32_t nnz_off_diag(std::int32_t row) const#

Number of non-zeros in unowned columns (off-diagonal block) on a given row.

Note

Can also be used on ghost rows

std::pair<std::span<const std::int32_t>, std::span<const std::int64_t>> graph() const#

Sparsity pattern graph after assembly. Uses local indices for the columns.

Note

Column global indices can be obtained from SparsityPattern::column_index_map()

Note

Includes ghost rows

Returns:

Adjacency list edges and offsets

std::span<const std::int32_t> off_diagonal_offsets() const#

Row-wise start of off-diagonals (unowned columns) for each row.

Note

Includes ghost rows

MPI_Comm comm() const#

Return MPI communicator.

template<typename T>
class SuperLUDistMatrix#
#include <superlu_dist.h>

SuperLU_DIST matrix interface.

Public Functions

SuperLUDistMatrix(const MatrixCSR<T> &A)#

Create SuperLU_DIST matrix operator.

Copies data from native CSR into SuperLU_DIST format.

Template Parameters:

T – Scalar type.

Parameters:

A – Matrix.

SuperLUDistMatrix(const SuperLUDistMatrix&) = delete#

Copy constructor (deleted).

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

Copy assignment (deleted).

MPI_Comm comm() const#

Get MPI communicator that matrix is defined on.

SuperLUDistStructs::SuperMatrix *supermatrix() const#

Get pointer to SuperLU_DIST SuperMatrix (non-const).

template<typename T>
class SuperLUDistSolver#
#include <superlu_dist.h>

SuperLU_DIST linear solver interface.

Public Functions

SuperLUDistSolver(std::shared_ptr<const SuperLUDistMatrix<T>> A)#

Create solver for a SuperLU_DIST matrix operator.

Solves linear system Au = b via LU decomposition.

The SuperLU_DIST solver has options set to upstream defaults, except PrintStat (verbose solver output) set to NO.

Template Parameters:

T – Scalar type.

Parameters:

A – Assembled left-hand side matrix.

SuperLUDistSolver(const SuperLUDistSolver&) = delete#

Copy constructor.

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

Copy assignment.

void set_option(std::string name, std::string value)#

Set solver option name to value.

See SuperLU_DIST User’s Guide for option names and values.

Parameters:

  • name – Option name.

  • value – Option value.

void set_options(SuperLUDistStructs::superlu_dist_options_t options)#

Set all solver options (native struct).

See SuperLU_DIST User’s Guide for option names and values.

Callers must complete the forward declared struct, e.g.:

#include <superlu_defs.h>
struct dolfinx::la::SuperLUDistStructs::superlu_dist_options_t
 : public ::superlu_dist_options_t
{
};

SuperLUDistStructs::superlu_dist_options_t options;
set_default_options_dist(&options);
options.PrintStat = YES;
// Setup SuperLUDistSolver
solver.set_options(options);
Parameters:

options – SuperLU_DIST option struct.

void set_A(std::shared_ptr<const SuperLUDistMatrix<T>> A)#

Set assembled left-hand side matrix A.

For advanced use with SuperLU_DIST option Factor allowing use of previously computed permutations when solving with new matrix A.

Parameters:

A – Assembled left-hand side matrix.

int solve(const Vector<T> &b, Vector<T> &u) const#

Solve linear system Au = b.

Note

Vectors must have size and parallel layout compatible with A.

Note

The caller must check the return code for success (== 0).

Note

The caller must u.scatter_forward() after the solve.

Note

The values of A are modified in-place during the solve.

Note

To solve with successive right-hand sides the caller must solver.set_options("Factor", "FACTORED") after the first solve.

Parameters:

  • b – Right-hand side vector.

  • u – Solution vector, overwritten during solve.

Returns:

SuperLU_DIST info integer.

class SuperLUDistStructs#
#include <superlu_dist.h>

Forward declare structs to avoid exposing SuperLU_DIST headers.

struct SuperMatrix : public SuperMatrix#
struct vec_int_t#

Struct holding vector of type int_t.

Public Members

std::vector<int_t> vec#

vector

struct gridinfo_t : public gridinfo_t#
struct superlu_dist_options_t : public superlu_dist_options_t#
struct dScalePermstruct_t : public dScalePermstruct_t#
struct sScalePermstruct_t : public sScalePermstruct_t#
struct zScalePermstruct_t : public zScalePermstruct_t#
struct dLUstruct_t : public dLUstruct_t#
struct sLUstruct_t : public sLUstruct_t#
struct zLUstruct_t : public zLUstruct_t#
struct dSOLVEstruct_t : public dSOLVEstruct_t#
struct sSOLVEstruct_t : public sSOLVEstruct_t#
struct zSOLVEstruct_t : public zSOLVEstruct_t#
struct SuperMatrixDeleter#
#include <superlu_dist.h>

Call library cleanup and delete pointer. For use with std::unique_ptr holding SuperMatrix.

Public Functions

void operator()(SuperLUDistStructs::SuperMatrix *A) const noexcept#

Deletion on SuperMatrix.

Parameters:

A –

template<typename T, typename Container = std::vector<T>, typename ScatterContainer = std::vector<std::int32_t>>
class Vector#
#include <Vector.h>

A vector that can be distributed across processes.

Template Parameters:

  • T – Scalar type of the vector.

  • Container – Data container type. This is typically std::vector<T> on CPUs, and thrust::device_vector<T> on GPUs.

  • ScatterContainer – Storage container type for the scatterer indices. This is typically std::vector<std::int32_t> on CPUs, and thrust::device_vector<std::int32_t> on GPUs.

Public Types

using container_type = Container#

Container type.

using value_type = container_type::value_type#

Scalar type.

Public Functions

inline Vector(std::shared_ptr<const common::IndexMap> map, int bs)#

Create a distributed vector.

Parameters:

  • map – Index map that describes the parallel layout of the data.

  • bs – Number of entries per index map ‘index’ (block size).

Vector(const Vector &x) = default#

Copy constructor.

Vector(Vector &&x) = default#

Move constructor.

template<typename T0, typename Container0, typename ScatterContainer0>
inline explicit Vector(const Vector<T0, Container0, ScatterContainer0> &x)#

Copy-convert vector, possibly using to different container types.

Examples of use include copying a Vector to a different value type, e.g. double to float, or copying a Vector from a CPU to a GPU.

Template Parameters:

Vec – Type of the Vector being copied.

Parameters:

x – Vector to copy.

Vector &operator=(Vector &&x) = default#

Move assignment operator.

inline void set(value_type v)#

Set all entries (including ghosts).

Deprecated:

Use std::ranges::fill(u.array(), v) instead.

Parameters:

v – [in] Value to set all entries to (on calling rank).

template<typename U, typename GetPtr>
inline void scatter_fwd_begin(U pack, GetPtr get_ptr)#

Begin scatter (send) of local data that is ghosted on other processes.

The user provides the function to pack to the send buffer. Typically usage would be a specialised function to pack data that resides on a GPU.

Note

Collective MPI operation

Template Parameters:

  • U – Pack function type.

  • GetPtr –

Parameters:

  • pack – Function that packs owned data into a send buffer.

  • get_ptr – Function that for a Container type returns the pointer to the underlying data.

inline void scatter_fwd_begin()#

Begin scatter (send) of local data that is ghosted on other processes (simplified CPU version).

Suitable for scatter operations on a CPU with std::vector storage. The send buffer is packed internally by a function that is suitable for use on a CPU.

Note

Collective MPI operation.

template<typename U>
inline void scatter_fwd_end(U unpack)#

End scatter (send) of local data values that are ghosted on other processes.

The user provides the function to unpack the receive buffer. Typically usage would be a specialised function to unpack data that resides on a GPU.

Note

Collective MPI operation.

Parameters:

unpack – Function to unpack the receive buffer into the ghost entries.

inline void scatter_fwd_end()#

End scatter (send) of local data values that are ghosted on other processes (simplified CPU version).

Suitable for scatter operations on a CPU with std::vector storage. The receive buffer is unpacked internally by a function that is suitable for use on a CPU.

Note

Collective MPI operation.

inline void scatter_fwd()#

Scatter (send) of local data values that are ghosted on other processes and update ghost entry values (simplified CPU version).

Suitable for scatter operations on a CPU with std::vector storage. The send buffer is packed and the receive buffer unpacked by a function that is suitable for use on a CPU.

Note

Collective MPI operation

template<typename U, typename GetPtr>
inline void scatter_rev_begin(U pack, GetPtr get_ptr)#

Start scatter (send) of ghost entry data to the owning process of an index.

The user provides the function to pack to the send buffer. Typically usage would be a specialised function to pack data that resides on a GPU.

Note

Collective MPI operation

Template Parameters:

  • U – Pack function type.

  • GetPtr –

Parameters:

  • pack – Function that packs ghost data into a send buffer.

  • get_ptr – Function that for a Container type returns the pointer to the underlying data.

inline void scatter_rev_begin()#

Start scatter (send) of ghost entry data to the owning process of an index (simplified CPU version).

Suitable for scatter operations on a CPU with std::vector storage. The send buffer is packed internally by a function that is suitable for use on a CPU.

Note

Collective MPI operation

template<typename U>
inline void scatter_rev_end(U unpack)#

End scatter of ghost data to owner and update owned entries.

For an owned entry, data from more than one process may be received. The received data can be summed or inserted into the owning entry by the unpack function.

Note

Collective MPI operation

template<class BinaryOperation>
inline void scatter_rev(BinaryOperation op)#

Scatter (send) of ghost data values to the owning process and assign/accumulate into the owned data entries (simplified CPU version).

For an owned entry, data from more than one process may be received. The received data can be summed or inserted into the owning entry. The this is controlled by the op function.

Note

Collective MPI operation

Parameters:

op – IndexMap operation (add or insert).

inline std::shared_ptr<const common::IndexMap> index_map() const#

Get IndexMap.

inline int bs() const#

Get block size.

inline container_type &array()#

Get the process-local part of the vector.

Owned entries appear first, followed by ghosted entries.

inline const container_type &array() const#

Get the process-local part of the vector (const version).

Owned entries appear first, followed by ghosted entries.

inline container_type &mutable_array()#

Get local part of the vector.

Deprecated:

Use ::array instead.

namespace impl#

Fetch the rows of B that correspond to the ghost columns of A.

For computing the product A*B, each rank needs the rows of B whose global indices match the ghost columns of A. Those ghost columns are owned by remote ranks, so we request those rows via neighbourhood communication.

Param A:

Matrix whose ghost column indices determine which rows of B are needed.

Param B:

Matrix whose rows are fetched.

Return:

A new MatrixCSR containing the fetched ghost rows of B, with an extended column IndexMap covering any new ghost columns introduced by those rows.

Functions

template<typename T>
std::tuple<std::shared_ptr<common::IndexMap>, std::vector<std::int32_t>, std::vector<std::int32_t>, std::vector<T>> fetch_ghost_rows(const dolfinx::la::MatrixCSR<T> &A, const dolfinx::la::MatrixCSR<T> &B)#

Fetch the rows of Matrix B which are referenced by the ghost columns of Matrix A.

Parameters:
Returns:

Tuple containing [new index map, rowptr, cols, values] for the received rows

template<typename T>
std::tuple<std::vector<std::int64_t>, std::vector<std::int32_t>, std::vector<std::int32_t>, std::vector<T>> matmul(const dolfinx::la::MatrixCSR<T> &A, const dolfinx::la::MatrixCSR<T> &B, std::shared_ptr<const common::IndexMap> new_col_map, std::span<const std::int32_t> ghost_row_ptr, std::span<const std::int32_t> ghost_cols, std::span<const T> ghost_vals)#

Compute the sparsity pattern and values of C = A*B in a single pass.

Uses a dense accumulator (one entry per possible output column) to simultaneously detect nonzero columns and accumulate A[i,j]*B[j,k], eliminating the separate value-fill loop and the per-entry lower_bound search that a two-pass approach requires.

Parameters:

  • A – Left matrix.

  • B – Right matrix (local rows only).

  • new_col_map – Extended column IndexMap for C, from fetch_ghost_rows.

  • ghost_row_ptr – CSR row pointer for the ghost rows of B.

  • ghost_cols – Local column indices (w.r.t. new_col_map) for ghost rows.

  • ghost_vals – Values for the ghost rows of B.

Returns:

Tuple (row_ptr, off_diag_offsets, cols, vals). off_diag_offsets[i] is the number of diagonal-block entries in row i, computed during the same sort step that establishes column order.

template<int BS0, int BS1, typename OP, typename U, typename V, typename W, typename X, typename Y>
void insert_csr(U &&data, const V &cols, const W &row_ptr, const X &x, const Y &xrows, const Y &xcols, OP op, typename Y::value_type num_rows)#

Incorporate data into a CSR matrix.

0 1 | 2 3 4 5 | 6 7

8 9 | 10 11 12 13 | 14 15

Note

In the case of block data, where BS0 or BS1 are greater than one, the layout of the input data is still the same. For example, the following can be inserted into the top-left corner with xrows={0,1} and xcols={0,1}, BS0=2, BS1=2 and x={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}.

Template Parameters:

  • BS0 – Row block size (of both matrix and data).

  • BS1 – Column block size (of both matrix and data).

  • OP – The operation (usually “set” or “add”).

Parameters:

  • data – [out] CSR matrix data.

  • cols – [in] CSR column indices.

  • row_ptr – [in] Pointer to the ith row in the CSR data.

  • x – [in] The m by n dense block of values (row-major) to add to the matrix.

  • xrows – [in] Row indices of x.

  • xcols – [in] Column indices of x.

  • op – [in] The operation (set or add),

  • num_rows – [in] Maximum row index that can be set. Used when debugging to check that rows beyond a permitted range are not being set.

template<int BS0, int BS1, typename OP, typename U, typename V, typename W, typename X, typename Y>
void insert_blocked_csr(U &&data, const V &cols, const W &row_ptr, const X &x, const Y &xrows, const Y &xcols, OP op, typename Y::value_type num_rows)#

Incorporate blocked data with given block sizes into a non-blocked MatrixCSR.

Note

Matrix block size (bs=1). Matrix sparsity must be correct to accept the data.

Note

See ::insert_csr for data layout.

Template Parameters:

  • BS0 – Row block size of data.

  • BS1 – Column block size of data.

  • OP – The operation (usually “set” or “add”).

Parameters:

  • data – [out] CSR matrix data.

  • cols – [in] CSR column indices.

  • row_ptr – [in] Pointer to the ith row in the CSR data.

  • x – [in] The m by n dense block of values (row-major) to add to the matrix.

  • xrows – [in] Row indices of x.

  • xcols – [in] Column indices of x.

  • op – [in] The operation (set or add).

  • num_rows – [in] Maximum row index that can be set. Used when debugging to check that rows beyond a permitted range are not being set.

template<typename OP, typename U, typename V, typename W, typename X, typename Y>
void insert_nonblocked_csr(U &&data, const V &cols, const W &row_ptr, const X &x, const Y &xrows, const Y &xcols, OP op, typename Y::value_type num_rows, int bs0, int bs1)#

Incorporate non-blocked data into a blocked matrix (data block size=1).

Note

Matrix sparsity must be correct to accept the data.

Note

See ::insert_csr for data layout.

Parameters:

  • data – [out] CSR matrix data.

  • cols – [in] CSR column indices.

  • row_ptr – [in] Pointer to the ith row in the CSR data.

  • x – [in] The m by n dense block of values (row-major) to add to the matrix.

  • xrows – [in] Rrow indices of x.

  • xcols – [in] Column indices of x.

  • op – [in] The operation (set or add).

  • num_rows – [in] Maximum row index that can be set. Used when debugging to check that rows beyond a permitted range are not being set.

  • bs0 – [in] Row block size of matrix.

  • bs1 – [in] Column block size of matrix.

template<typename T, int BS1>
void spmv(std::span<const T> values, std::span<const std::int64_t> row_begin, std::span<const std::int64_t> row_end, std::span<const std::int32_t> indices, std::span<const T> x, std::span<T> y, int bs0, int bs1)#

Sparse matrix-vector product implementation.

Template Parameters:

  • T –

  • BS1 –

Parameters:

  • values –

  • row_begin –

  • row_end –

  • indices –

  • x –

  • y –

  • bs0 –

  • bs1 –

template<typename T, int BS1>
void spmvT(std::span<const T> values, std::span<const std::int64_t> row_begin, std::span<const std::int64_t> row_end, std::span<const std::int32_t> indices, std::span<const T> x, std::span<T> y, int bs0, int bs1)#

Sparse matrix-vector transpose product implementation.

Computes y += A^T x where A is given in CSR format. The transpose is applied implicitly: for each nonzero A(i, indices[j]) the contribution A(i,j) * x[i] is scattered into y[indices[j]].

Note

y is accumulated into (not overwritten). Callers should zero y before the first call if a fresh result is required.

Note

The value block layout is row-major within each block: for block entry j the element at row-offset k0 and column-offset k1 is stored at values[j * bs0 * bs1 + k0 * bs1 + k1].

Template Parameters:

  • T – Scalar type of the matrix and vector entries.

  • BS1 – Compile-time column block size. Pass -1 to use the runtime value bs1 instead.

Parameters:

  • values – [in] Nonzero values of A, stored block-row-major. Length: nnz * bs0 * bs1.

  • row_begin – [in] Start positions in values/indices for each row of A. Length: number of rows of A.

  • row_end – [in] End positions in values/indices for each row of A. Length: number of rows of A.

  • indices – [in] Column indices of each nonzero block entry of A. Length: nnz.

  • x – [in] Input vector, indexed by the rows of A. Length: num_rows * bs0.

  • y – [inout] Output vector, indexed by the columns of A, accumulated into. Length: num_cols * bs1.

  • bs0 – [in] Row block size (runtime value).

  • bs1 – [in] Column block size (runtime value, used when BS1 == -1).

template<typename T, int BS0 = -1, int BS1 = -1>
std::tuple<std::vector<std::int32_t>, std::vector<std::int64_t>, std::vector<T>> local_transpose(const dolfinx::la::MatrixCSR<T> &A)#

Compute the local (diagonal-block) transpose of A.

Iterates over owned rows of A and the owned-column entries within each row (indices [row_ptr[i], off_diag_offset[i])), building the transpose CSR in one bucket-fill pass. The resulting column indices are original row indices (0-based local), so columns within every output row are already in ascending order.

Parameters:

A – MatrixCSR

Template Parameters:

  • T – Scalar type

  • BS0 – row block size, must match row block size of A, or use -1 for non-optimized

  • BS1 – col block size, must match col block size of A, or use -1 for non-optimized

Returns:

Tuple (cols, row_ptr, values) for the transposed diagonal block.

struct Sparsity#
#include <matmul.h>

Lightweight sparsity descriptor satisfying the SparsityImplementation concept required by the MatrixCSR constructor.

Holds non-owning views into externally managed CSR arrays together with the row and column IndexMaps that describe the parallel distribution. It is intended as a short-lived helper: to be constructed it from the output of impl::matmul and passed to the MatrixCSR constructor.

Note

All spans must remain valid for the lifetime of this object.

Public Functions

inline std::shared_ptr<const common::IndexMap> index_map(int dim) const#

Return the row (dim == 0) or column (dim == 1) IndexMap.

inline int block_size(int i) const#

Return the block size.

inline std::pair<std::span<const std::int32_t>, std::span<const std::int64_t>> graph() const#

Return the CSR graph as (column_indices, row_pointers).

inline std::span<const std::int32_t> off_diagonal_offsets() const#

Return the per-row diagonal block sizes (see _off_diag).

namespace petsc#

PETSc linear algebra functions.

Functions

void error(PetscErrorCode error_code, const std::string &filename, const std::string &petsc_function)#

Print error message for PETSc calls that return an error.

std::vector<Vec> create_vectors(MPI_Comm comm, const std::vector<std::span<const PetscScalar>> &x)#

Create PETsc vectors from the local data. The data is copied into the PETSc vectors and is not shared.

Note

Caller is responsible for destroying the returned object

Parameters:

  • comm – [in] The MPI communicator

  • x – [in] The vector data owned by the calling rank. All components must have the same length.

Returns:

Array of PETSc vectors

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

Create a ghosted PETSc Vec

Note

Caller is responsible for destroying the returned object

Parameters:

  • map – [in] The index map describing the parallel layout (by block)

  • bs – [in] The block size

Returns:

A PETSc Vec

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

Create a ghosted PETSc Vec from a local range and ghost indices

Note

Caller is responsible for freeing the returned object

Parameters:

  • comm – [in] The MPI communicator

  • range – [in] The local ownership range (by blocks)

  • ghosts – [in] Ghost blocks

  • bs – [in] The block size. The total number of local entries is bs * (range[1] - range[0]).

Returns:

A PETSc Vec

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

Create a PETSc Vec that wraps the data in an array

Note

The array x must be kept alive to use the PETSc Vec object

Note

The caller should call VecDestroy to free the return PETSc vector

Parameters:

  • map – [in] The index map that describes the parallel layout of the distributed vector (by block)

  • bs – [in] Block size

  • x – [in] The local part of the vector, including ghost entries

Returns:

A PETSc Vec object that shares the data in x

template<class V>
Vec create_vector_wrap(const la::Vector<V> &x)#

Create a PETSc Vec that wraps the data in an array

Parameters:

x – [in] The vector to be wrapped

Returns:

A PETSc Vec object that shares the data in x

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

Compute PETSc IndexSets (IS) for a stack of index maps.

If map[0] = {0, 1, 2, 3, 4, 5, 6} and map[1] = {0, 1, 2, 4} (in local indices) then IS[0] = {0, 1, 2, 3, 4, 5, 6} and IS[1] = {7, 8, 9, 10}.

Todo:

This function could take just the local sizes.

Note

The caller is responsible for destruction of each IS.

Parameters:

maps – [in] Vector of IndexMaps and corresponding block sizes

Returns:

Vector of PETSc Index Sets, created on PETSC_COMM_SELF

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

Copy blocks from Vec into local arrays.

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

Scatter local vectors to Vec.

Mat create_matrix(MPI_Comm comm, const SparsityPattern &sp, std::optional<std::string> type = std::nullopt)#

Create a PETSc Mat. Caller is responsible for destroying the returned object.

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

Create PETSc MatNullSpace. Caller is responsible for destruction returned object.

Parameters:

  • comm – [in] The MPI communicator

  • basis – [in] The nullspace basis vectors

Returns:

A PETSc nullspace object

class Vector#
#include <petsc.h>

A simple wrapper for a PETSc vector pointer (Vec). Its main purpose is to assist with memory/lifetime management of PETSc Vec objects.

Access the underlying PETSc Vec pointer using the function Vector::vec() and use the full PETSc interface.

Public Functions

Vector(const common::IndexMap &map, int bs)#

Create a vector

Note

Collective

Parameters:

  • map – [in] Index map describing the parallel layout

  • bs – [in] the block size

Vector(Vector &&x) noexcept#

Move constructor.

Vector(Vec x, bool inc_ref_count)#

Create holder of a PETSc Vec object/pointer. The Vec x object should already be created. If inc_ref_count is true, the reference counter of the Vec object will be increased. The Vec reference count will always be decreased upon destruction of the PETScVector.

Note

Collective

Parameters:

  • x – [in] The PETSc Vec

  • inc_ref_count – [in] True if the reference count of x should be incremented

virtual ~Vector()#

Destructor.

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

Move Assignment operator.

Vector copy() const#

Create a copy of the vector

Note

Collective

std::int64_t size() const#

Return global size of the vector.

std::int32_t local_size() const#

Return local size of vector (belonging to the call rank).

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

Return ownership range for calling rank.

MPI_Comm comm() const#

Return MPI communicator.

void set_options_prefix(const std::string &options_prefix)#

Sets the prefix used by PETSc when searching the options database.

std::string get_options_prefix() const#

Returns the prefix used by PETSc when searching the options database

void set_from_options()#

Call PETSc function VecSetFromOptions on the underlying Vec object.

Vec vec() const#

Return pointer to PETSc Vec object.

class Operator#
#include <petsc.h>

This class is a base class for matrices that can be used in petsc::KrylovSolver.

Subclassed by Matrix

Public Functions

Operator(Mat A, bool inc_ref_count)#

Constructor.

Operator(Operator &&A) noexcept#

Move constructor.

virtual ~Operator()#

Destructor.

Operator &operator=(const Operator &A) = delete#

Assignment operator (deleted).

Operator &operator=(Operator &&A) noexcept#

Move assignment operator.

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

Return number of rows and columns (num_rows, num_cols). PETSc returns -1 if size has not been set.

Vec create_vector(std::size_t dim) const#

Initialize vector to be compatible with the matrix-vector product y = Ax. In the parallel case, size and layout are both important.

Parameters:

dim – [in] The dimension (axis): dim = 0 –> z = y, dim = 1 –> z = x

Mat mat() const#

Return PETSc Mat pointer.

class Matrix : public Operator#
#include <petsc.h>

It is a simple wrapper for a PETSc matrix pointer (Mat). Its main purpose is to assist memory management of PETSc Mat objects.

For advanced usage, access the PETSc Mat pointer using the function mat() and use the standard PETSc interface.

Public Types

enum class AssemblyType : std::int8_t#

Assembly type FINAL - corresponds to PETSc MAT_FINAL_ASSEMBLY FLUSH - corresponds to PETSc MAT_FLUSH_ASSEMBLY

Values:

enumerator FINAL#
enumerator FLUSH#

Public Functions

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

Create holder for a PETSc Mat object from a sparsity pattern.

Matrix(Mat A, bool inc_ref_count)#

Create holder of a PETSc Mat object/pointer. The Mat A object should already be created. If inc_ref_count is true, the reference counter of the Mat will be increased. The Mat reference count will always be decreased upon destruction of the petsc::Matrix.

Matrix(Matrix &&A) = default#

Move constructor (falls through to base class move constructor).

~Matrix() = default#

Destructor.

Matrix &operator=(const Matrix &A) = delete#

Assignment operator (deleted).

Matrix &operator=(Matrix &&A) = default#

Move assignment operator.

void apply(AssemblyType type)#

Finalize assembly of tensor. The following values are recognized for the mode parameter:

Parameters:

type – FINAL - corresponds to PETSc MatAssemblyBegin+End(MAT_FINAL_ASSEMBLY) FLUSH - corresponds to PETSc MatAssemblyBegin+End(MAT_FLUSH_ASSEMBLY)

double norm(Norm norm_type) const#

Return norm of matrix.

void set_options_prefix(const std::string &options_prefix)#

Sets the prefix used by PETSc when searching the options database

std::string get_options_prefix() const#

Returns the prefix used by PETSc when searching the options database

void set_from_options()#

Call PETSc function MatSetFromOptions on the PETSc Mat object.

Public Static Functions

static inline auto set_fn(Mat A, InsertMode mode)#

Return a function with an interface for adding or inserting values into the matrix A (calls MatSetValuesLocal)

Parameters:

  • A – [in] The matrix to set values in

  • mode – [in] The PETSc insert mode (ADD_VALUES, INSERT_VALUES, …)

static inline auto set_block_fn(Mat A, InsertMode mode)#

Return a function with an interface for adding or inserting values into the matrix A using blocked indices (calls MatSetValuesBlockedLocal)

Parameters:

  • A – [in] The matrix to set values in

  • mode – [in] The PETSc insert mode (ADD_VALUES, INSERT_VALUES, …)

static inline auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)#

Return a function with an interface for adding or inserting blocked values to the matrix A using non-blocked insertion (calls MatSetValuesLocal). Internally it expands the blocked indices into non-blocked arrays.

Parameters:

  • A – [in] The matrix to set values in

  • bs0 – [in] Block size for the matrix rows

  • bs1 – [in] Block size for the matrix columns

  • mode – [in] The PETSc insert mode (ADD_VALUES, INSERT_VALUES, …)

class KrylovSolver#
#include <petsc.h>

This class implements Krylov methods for linear systems of the form Ax = b. It is a wrapper for the Krylov solvers of PETSc.

Public Functions

explicit KrylovSolver(MPI_Comm comm)#

Create Krylov solver for a particular method and named preconditioner

KrylovSolver(KSP ksp, bool inc_ref_count)#

Create solver wrapper of a PETSc KSP object

Parameters:

  • ksp – [in] The PETSc KSP object. It should already have been created

  • inc_ref_count – [in] Increment the reference count on ksp if true

KrylovSolver(KrylovSolver &&solver) noexcept#

Move constructor.

~KrylovSolver()#

Destructor.

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

Move assignment.

void set_operator(const Mat A)#

Set operator (Mat).

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

Set operator and preconditioner matrix (Mat).

int solve(Vec x, const Vec b, bool transpose = false) const#

Solve linear system Ax = b and return number of iterations (A^t x = b if transpose is true)

void set_options_prefix(const std::string &options_prefix)#

Sets the prefix used by PETSc when searching the PETSc options database

std::string get_options_prefix() const#

Returns the prefix used by PETSc when searching the PETSc options database

void set_from_options() const#

Set options from PETSc options database.

KSP ksp() const#

Return PETSc KSP pointer.

namespace options#

These class provides static functions that permit users to set and retrieve PETSc options via the PETSc option/parameter system. The option must not be prefixed by ‘-’, e.g. 

la::petsc::options::set("mat_mumps_icntl_14", 40);

Functions

void set(std::string option)#

Set PETSc option that takes no value.

template<typename T>
void set(std::string option, const T &value)#

Generic function for setting PETSc option.

void clear(std::string option)#

Clear a PETSc option.

void clear()#

Clear PETSc global options database.

On this page
Show Source