GitHub

Mat

PETSc matrices (Mat) provide sparse and dense matrix storage with efficient parallel operations. They are essential for discretizing PDEs and setting up linear/nonlinear systems.

Overview

PETSc matrices support:

  • Sparse formats: AIJ (CSR), BAIJ (block CSR), and more
  • Dense format: For small matrices or dense operations
  • Parallel distribution: Row-based distribution across MPI processes
  • Matrix-free operations: Via MatShell for custom operators

PetscMat is the one constructor: it replaces v0.4's MatSeqAIJ, MatSeqDense, MatCreateSeqAIJ, MatSeqAIJWithArrays and MatAIJ (naming conventions, §6). The old spellings still work in v0.5 and warn once.

Creating Matrices

Sparse Matrices (AIJ/CSR Format)

# Create sparse matrix with estimated non-zeros per row
A = PetscMat(petsclib, num_rows, num_cols, nnz_per_row)

# From Julia SparseMatrixCSC
using SparseArrays
S = sprand(100, 100, 0.1)
A = PetscMat(petsclib, MPI.COMM_SELF, S)

# With varying non-zeros per row
nnz = petsclib.PetscInt[5, 3, 4]  # One value per row
A = PetscMat(petsclib, num_rows, num_cols, nnz)

Dense Matrices

# Wrap a Julia matrix (no copy)
julia_mat = rand(10, 10)
A = PetscMat(petsclib, julia_mat)

From DM Objects

# Create matrix with sparsity pattern from DM
A = LibPETSc.DMCreateMatrix(petsclib, dm)

Matrix Shell (Matrix-Free)

# Create a shell matrix with custom mult operation, y = mult_function(y, x)
A = PETSc.MatShell(petsclib, mult_function, MPI.COMM_SELF, local_rows, local_cols)

Setting Values

# Set individual element. `setindex!` is 1-based, and converts for you
A[i, j] = value

# `set_values!` is the bulk route, and it keeps PETSc's 0-based indices: the
# `_0b` in the parameter names says so (naming conventions, §12.1)
PETSc.set_values!(A, rows_0b, cols_0b, values, PETSc.INSERT_VALUES)

# The same call with `MatStencil` rows and columns, for stencil-based assembly
PETSc.set_values!(A, stencils_0b, stencils_0b, values, PETSc.INSERT_VALUES)

Assembly

Matrices must be assembled after setting values:

# Set all values first
A[1, 1] = 2.0
A[1, 2] = -1.0
# ...

# Then assemble
PETSc.assemble!(A)

Common Operations

size(A)                    # Get (rows, cols)
PETSc.ownership_range(A)   # The 1-based rows owned by this process
PETSc.setup!(A)            # Complete matrix setup
PETSc.destroy!(A)          # Release it

ownership_range(A) is 1-based and takes no second argument. v0.4's ownershiprange(A, false) still works and warns; it is a MethodError in v0.6 (naming conventions, §12.1).

Functions

PETSc.LibPETSc.PetscMat — Method
M::PetscMat = PetscMat(petsclib, S::SparseMatrixCSC; with_arrays = false)
M::PetscMat = PetscMat(petsclib, comm, S::SparseMatrixCSC; with_arrays = false)

Creates a PetscMat object from a Julia SparseMatrixCSC S in sequential AIJ format. comm defaults to MPI.COMM_SELF.

By default the entries are copied into storage PETSc allocates. With with_arrays = true the CSR arrays converted from S are handed to PETSc and borrowed rather than copied, which is what v0.4's MatSeqAIJWithArrays did; the arrays are kept alive for as long as the matrix.

Replaces v0.4's MatCreateSeqAIJ: construction goes through the type (docs/src/man/naming.md §5.1). The unrelated MatSeqAIJ, which allocated from sizes, is now PetscMat(petsclib, m, n, nnz).

External Links

source
PETSc.LibPETSc.PetscMat — Method
B = PetscMat(petsclib, rowptr, colval, nzval; comm = MPI.COMM_SELF, ncols = …)

Create a PETSc SeqAIJ matrix directly on the CSR arrays rowptr, colval and nzval, which PETSc borrows rather than copies.

rowptr and colval are 0-based, PETSc's own base for bulk index arrays (docs/src/man/naming.md §12.1). The number of rows is length(rowptr) - 1; ncols defaults to one past the largest column index.

The matrix keeps the arrays alive for as long as it exists, including when a solver still holds it after destroy! on this handle.

Replaces v0.4's MatSeqAIJWithArrays, which took a SparseMatrixCSC and so could not be told apart from MatCreateSeqAIJ (docs/src/man/naming.md §6).

External Links

source
PETSc.LibPETSc.PetscMat — Method
mat = PetscMat(petsclib, num_rows, num_cols, nonzeros)

Create a PETSc serial sparse array using AIJ format (also known as a compressed sparse row or CSR format) of size num_rows X num_cols with nonzeros per row

If nonzeros is an Integer the same number of non-zeros will be used for each row, if nonzeros is a Vector{PetscInt} then one value must be specified for each row.

Memory allocation is handled by PETSc and garbage collection can be used.

External Links

source
PETSc.LibPETSc.PetscMat — Method
mat = PetscMat(petsclib, num_rows, num_cols; type = :seqaij)

An empty PETSc matrix of the given size.

type is the PETSc implementation name as a Symbol (docs/src/man/naming.md §3.1): :seqaij (the default) preallocates nothing, :seqdense (spelled :dense as well) allocates the dense storage. Flavour is a keyword rather than a type parameter because it only affects construction, and every variant returns a PetscMat (§9).

External Links

source
PETSc.LibPETSc.PetscMat — Method
mat = PetscMat(petsclib, A::Matrix{PetscScalar})

PETSc dense array. This wraps a Julia Matrix{PetscScalar} object.

Replaces v0.4's MatSeqDense: construction goes through the type (docs/src/man/naming.md §5.1).

External Links

source
PETSc.MatPtr — Type
MatPtr(petsclib, ptr::CMat, own::Bool)

Container type for a PETSc Mat that is just a raw pointer.

If own is true a finalizer is set on the matrix, but only on a serial communicator, since MatDestroy is collective and a GC finalizer runs at an arbitrary point. If own is false the handle belongs to PETSc and destroy! is a no-op, leaving the wrapper usable.

source
PETSc.MatShell — Type
MatShell(
    petsclib::PetscLib,
    obj::OType,
    comm::MPI.Comm,
    local_rows,
    local_cols,
    global_rows = LibPETSc.PETSC_DECIDE,
    global_cols = LibPETSc.PETSC_DECIDE,
)

Create a global_rows X global_cols PETSc shell matrix object wrapping obj with local size local_rows X local_cols.

The obj will be registered as an MATOP_MULT function and if if obj is a Function, then the multiply action obj(y,x); otherwise it calls mul!(y, obj, x).

if comm == MPI.COMM_SELF then the garbage connector can finalize the object, otherwise the user is responsible for calling destroy!.

Callback

The closure is kept with the PETSc object, not with the wrapper passed here, so a borrowed handle is fine to pass. It succeeds by returning and fails by throwing; its return value is ignored, except that a nonzero Integer still fails the call and warns until v0.6. An exception it throws comes out of the solve!, step! or setup! that ran it.

External Links

source
Base.fill! — Method
fill!(A::AbstractPetscMat, 0)

Set every stored entry of A to zero and return A. The nonzero pattern stays, so the matrix can be refilled without a new allocation. Only zero is accepted: PETSc has no operation that sets every entry to another value, and any other x throws an ArgumentError.

External Links

source
LinearAlgebra.mul! — Method
mul!(y::PetscVec{PetscLib}, M::AbstractPetscMat{PetscLib}, x::PetscVec{PetscLib})

Computes y = M*x

source
PETSc.csr_from_csc — Method
rowptr, colval, nzval = csr_from_csc(petsclib, A::SparseMatrixCSC)

The CSR triple PETSc wants, built from Julia's CSC storage. rowptr and colval are 0-based, as MatCreateSeqAIJWithArrays expects.

source
PETSc.destroy! — Method
destroy!(m::AbstractPetscMat)

Destroy a Mat (matrix) object and release associated resources.

This function is typically called automatically via finalizers when the object is garbage collected, but can be called explicitly to free resources immediately. Does nothing on a matrix that only borrows its handle: see owns.

External Links

source
PETSc.diagonal! — Method
diagonal!(d::AbstractPetscVec, A::AbstractPetscMat)

Write the diagonal of A into d and return d. d needs the row layout of A, as the left vector from MatCreateVecs has. LinearAlgebra.diag is not extended, because it returns a new vector.

External Links

source
PETSc.mat_seqaij_with_arrays — Method
mat_seqaij_with_arrays(petsclib, comm, A::SparseMatrixCSC)

The v0.4 MatSeqAIJWithArrays body, kept for its deprecation shim: converts A to CSR and hands the arrays to the PetscMat constructor.

source
PETSc.ownership_range — Method
ownership_range(mat::AbstractPetscMat)

The range of row indices owned by this processor, assuming that the mat is laid out with the first n1 rows on the first processor, next n2 rows on the second, etc. For certain parallel layouts this range may not be well defined.

The range is 1-based, always: an index into Julia data is 1-based (docs/src/man/naming.md §12.1). v0.4 took base_one::Bool positionally and made the convention a runtime choice; ownership_range(A, false) warns in v0.5 and is a MethodError in v0.6.

Note

unlike the C function, the range returned is inclusive (idx_first:idx_last)

External Links

source
PETSc.set_option! — Method
set_option!(A::AbstractPetscMat, option::LibPETSc.MatOption, flag::Bool)

Turn the matrix option option on or off and return A, e.g. set_option!(A, LibPETSc.MAT_NEW_NONZERO_ALLOCATION_ERR, false) to allow insertions outside the preallocated pattern.

External Links

source
PETSc.set_type! — Method
set_type!(A::AbstractPetscMat, type::Symbol)

Set the matrix implementation, for example :seqaij or :dense.

External Links

source
PETSc.set_values! — Method
set_values!(
    M::AbstractMat{PetscLib},
    rows_0b::Vector{PetscInt},
    cols_0b::Vector{PetscInt},
    rowvals::Array{PetscScalar},
    insertmode::InsertMode = INSERT_VALUES;
    num_rows = length(rows_0b),
    num_cols = length(cols_0b)
)

Set values of the matrix M with base-0 row and column indices rows_0b and cols_0b, inserting the values rowvals.

The _0b suffix says what the base is: a bulk index array handed to C keeps PETSc's base rather than being rebuilt on a hot path (docs/src/man/naming.md §12.1). A[i, j] = v is the 1-based route.

If the keyword arguments num_rows or num_cols is specified then only the first num_rows * num_cols values of rowvals will be used.

External Links

source
PETSc.set_values! — Method
set_values!(
    M::AbstractPetscMat{PetscLib},
    rows_0b::Vector{MatStencil},
    cols_0b::Vector{MatStencil},
    rowvals::Array{PetscScalar},
    insertmode::InsertMode = INSERT_VALUES;
    num_rows = length(rows_0b),
    num_cols = length(cols_0b)
)

Set values of the matrix M with base-0 row and column indices rows_0b and cols_0b, inserting the values rowvals.

The _0b suffix says what the base is: a bulk index array handed to C keeps PETSc's base rather than being rebuilt on a hot path (docs/src/man/naming.md §12.1). A[i, j] = v is the 1-based route.

If the keyword arguments num_rows or num_cols is specified then only the first num_rows * num_cols values of rowvals will be used.

External Links

source
PETSc.type_name — Method
type_name(A::AbstractPetscMat)

The name PETSc knows this matrix's implementation by, as a Symbol (:seqaij, :mpiaij, …), or nothing when the matrix has no type yet (docs/src/man/naming.md §3.1). v0.4 answered with a String, using the sentinel "(not set)"; that is a break with no shim (§16).

External Links

source
PETSc.zero_rows! — Method
zero_rows!(A::AbstractPetscMat, rows_0b, diag = 1; x = nothing, b = nothing)

Zero the rows rows_0b of the assembled matrix A, put diag on their diagonal entries and return A. With x and b given, also set b[i] = diag * x[i] for each zeroed row i, so that a solve keeps the values of x there: the usual way to impose a Dirichlet condition.

rows_0b are 0-based global row numbers, as PETSc takes them (docs/src/man/naming.md §12.1). Collective: every process calls it, each with its own rows, possibly none. zero_rows_local! takes local numbers instead.

External Links

source
PETSc.zero_rows_local! — Method
zero_rows_local!(A::AbstractPetscMat, rows_0b, diag = 1; x = nothing, b = nothing)

zero_rows! with 0-based local row numbers, translated through the local-to-global mapping of A. A matrix from DMCreateMatrix has that mapping; one built without it throws a PetscError.

External Links

source