Linear algebra (dolfinx::la)

namespace la

Linear algebra interface.

Interface to linear algebra data structures and solvers.

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

Norm types.

Values:

enumerator l1
enumerator l2
enumerator linf
enumerator frobenius

Functions

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

Compute the inner product of two vectors. The two vectors must have the same parallel layout

Note

Collective MPI operation

Parameters:

  • a – A vector

  • b – A vector

Returns:

Returns a^{H} b (a^{T} b if a and b are real)

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 – A vector

  • type – Norm type

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

Orthonormalize a set of vectors

Parameters:

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

Template Parameters:

V – dolfinx::la::Vector

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

Test if basis is orthonormal.

Returns true if ||x_i - x_j|| - delta_{ij} < eps fro 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.

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

Distributed sparse matrix.

The matrix storage format is compressed sparse row. The matrix is partitioned row-wise across MPI ranks.

Warning

Experimental storage of a matrix in CSR format which can be assembled into using the usual DOLFINx assembly routines. Matrix internal data can be accessed for interfacing with other code.

Template Parameters:

  • Scalar – Scalar type of matrix entries

  • Container – Sequence container type to store 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

MatrixCSR(const SparsityPattern &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

inline void set(value_type x)

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

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.

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.

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().

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.

double squared_norm() const

Compute the Frobenius norm squared across all processes.

Note

MPI Collective

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

Index maps for the row and column space.

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

Returns:

Row (0) or column (1) index maps

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 manually if not required.

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

Block size

Returns:

block sizes for rows and columns

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, const std::array<std::shared_ptr<const common::IndexMap>, 2> &maps, const 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.

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.

common::IndexMap column_index_map() const

Builds the index map for columns after assembly of the sparsity pattern.

Todo:

Should this be compted and stored when finalising the SparsityPattern?

Returns:

Map for all 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, typename Container = std::vector<T>>
class Vector
#include <Vector.h>

Distributed vector

Template Parameters:

  • T – Scalar type

  • Container – data container type

Public Types

using value_type = T

Scalar type.

using container_type = Container

Container type.

Public Functions

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

Create a distributed vector

Parameters:

  • map – IndexMap for parallel distribution of the data

  • bs – Block size

inline Vector(const Vector &x)

Copy constructor.

inline Vector(Vector &&x)

Move constructor.

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

Move Assignment operator.

inline void set(value_type v)

Set all entries (including ghosts)

Parameters:

v – [in] The value to set all entries to (on calling rank)

inline void scatter_fwd_begin()

Begin scatter of local data from owner to ghosts on other ranks

Note

Collective MPI operation

inline void scatter_fwd_end()

End scatter of local data from owner to ghosts on other ranks

Note

Collective MPI operation

inline void scatter_fwd()

Scatter local data to ghost positions on other ranks

Note

Collective MPI operation

inline void scatter_rev_begin()

Start scatter of ghost data to owner

Note

Collective MPI operation

template<class BinaryOperation>
inline void scatter_rev_end(BinaryOperation op)

End scatter of ghost data to owner. This process may receive data from more than one process, and the received data can be summed or inserted into the local portion of the vector.

Note

Collective MPI operation

Parameters:

op – The operation to perform when adding/setting received values (add or insert)

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

Scatter ghost data to owner. This process may receive data from more than one process, and the received data can be summed or inserted into the local portion of the vector.

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 std::span<const value_type> array() const

Get local part of the vector (const version)

inline std::span<value_type> mutable_array()

Get local part of the vector.

namespace impl

Functions

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] The CSR matrix data

  • cols – [in] The CSR column indices

  • row_ptr – [in] The 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] The row indices of x

  • xcols – [in] The column indices of x

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

  • num_rows – [in] The 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] The CSR matrix data

  • cols – [in] The CSR column indices

  • row_ptr – [in] The 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] The row indices of x

  • xcols – [in] The column indices of x

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

  • num_rows – [in] The 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] The CSR matrix data

  • cols – [in] The CSR column indices

  • row_ptr – [in] The 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] The row indices of x

  • xcols – [in] The column indices of x

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

  • num_rows – [in] The 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