Linear Operators

This section describes the different catagories of linear operators, including the API and the requirements for defining new linear operators. Note that whenever a linear solver is used, the Linear Solver interface should be used rather than the Linear Operator interface described here.

Type Hierarchy

When creating a new linear operator, it is very important to make it inherit from the proper abstract type.

LinearSolvers.AbstractLO — Type.

Abstract supertype of all linear operators used for A when solving Ax = b.

The purpose of this type is to provide a consistent interface for different types of linear operators. This type really combines two notions: what the type of the linear operator is, and how it should be solved.

Any implementation of this type should subtype the appropriate catagory of: AbstractDenseLO, AbstractSparseDirectLO, AbstractPetscMatLO, AbstractPetscMatFreeLO

Note that matrix-explicit implementations can often write a single function for all these cases if using an matrix interface functions that are defined for all the matrix types. See MatExplicitLO

Required Fields

lo_inner: another AbstractPC. Can be one of DenseLO, SparseDirectLO, PetscMatLO, or PetscMatFreeLO, or any other user defined linear operator.

LinearSolvers.AbstractDenseLO — Type.

Linear operator type for Dense matrices. This is generally used only for debugging.

LinearSolvers.AbstractSparseDirectLO — Type.

Linear operator type for SparseMatrixCSC matrices, which use a direct solver.

LinearSolvers.AbstractPetscMatLO — Type.

Linear operator type for Petsc matrix-explicit.

LinearSolvers.AbstractPetscMatFreeLO — Type.

Linear operator type for Petsc matrix-free.

Several typealiases are also available:

LinearSolvers.MatExplicitLO — Constant.

Useful union for all the matrix-explicit linear operator types. Because matrices have a small set of common interface functions, it is often possible to write a single function that works on all the different types of matrices.

LinearSolvers.PetscLO — Constant.

Union of Petsc linear operator types

LinearSolvers.DirectLO — Constant.

Union of linear operators that do direct solves

API

Every linear operator supports the following functions. When defining a new linear operator, some of the functions must be extended with new methods, while others will be created automatically based on the supertype of the linear operator.

LinearSolvers.calcLinearOperator — Function.

This function calculates the linear operator. Every implementation of AbstractLO should extend this function with a new method. For matrix-free operators, this function must exist but need not perform any actions.

For matrix-explicit implementations, this function should be called on lo_inner first and modifications to the Jacobian made subsequently.

Inputs

lo: the AbstractLO implementation (fields may be updated)
mesh
sbp
eqn
opts
ctx_residual: the ctx required by physicsRhs
t: current time

Implementation Notes: For matrix-free operation, this function sets the Petsc ctx for the PetscMat, which contains a reference to the mesh, sbp, eqn, opts arguments. This could lead to unexpected behavior if those arguments are modified and this function is not called again before the next solve.

Calculates the linear operator. Use this function only if you want to calculate the linear operator and not the preconditioner. Prefer calcPCandLO, which avoids calculating the matrix twice if the preconditioner and linear operator share the same matrix

Inputs

ls: StandardLinearSolver
mesh
sbp
eqn
opts
ctx_residual: the ctx required by physicsRhs like functions
t: current time

Keyword Arguments

start_comm: start parallel communication (if required by the lo), default false. This means the user is generally required to make sure parallel communication is started before calling this function.

LinearSolvers.applyLinearOperator — Function.

Applies the linear operator, ie. , Ax = b

Matrix-explicit implementations AbstractLO do not have to implement this function, though matrix-free implementations must extend it with a new method.

Inputs

lo: the AbstractLO implementation.
mesh
sbp
eqn
opts
ctx_residual: the ctx for physicsRhs or the another right hand side function built on top of it
t: current time
x: an AbstractVector (although never a PetscVec)

Inputs/Outputs

b: vector updated with results (do not overwrite)

LinearSolvers.applyLinearOperatorTranspose — Function.

Applies the transpose of the linear operator, ie. A.'*x = b Similar to applyLinearOperator, see that function for details.

Note that not every method supports this. In particular, Petsc matrix-free LinearOperators don't currently expose this (although they could with enough reverse-mode)

LinearSolvers.getBaseLO — Function.

Similar to getBasePC except it gets the underlying linear operator, ie. one of DenseLO, SparseDirectLO, PetscMatLO or PetscMatFreeLO.

For matrix-explicit methods, this is a good way of getting the underlying linear operator object, which contains the matrix in the A field (for all matrix-explicit linear operators).

Inputs

lo: an AbstractLO

LinearSolvers.getBaseObject — Function.

This function calls either getBasePC or getBaseLO depending on the type of its argument.

Inputs

lo: either an AbstractLO or AbstractPC object

Outputs

the base PC or LO object

Utils.free — Method.

This function frees any memory belonging to external libraries. Users must call this function when they are finished with an AbstractLO object.

Users do not have to define this function for their AbstractLO types.

Inputs

lo: the AbstractLO object

Concrete LO Types

LinearSolvers.DenseLO — Type.

Dense array linear operator. Serial only.

Public Fields

A: an Array{Float64, 1}

LinearSolvers.DenseLO — Method.

Outer constructor for DenseLO

Inputs

pc: a PCNone
mesh
sbp
eqn
opts

LinearSolvers.SparseDirectLO — Type.

LinearOperator type to be used with sparse direct solve. Note that the sparsity pattern of the matrix A must be constant.

Fields

A: the matrix

LinearSolvers.SparseDirectLO — Method.

Outer constructor for SparseDirectLO

Inputs

pc: a PCNone
mesh
sbp
eqn
opts

LinearSolvers.PetscMatLO — Type.

Petsc matrix-explicit linear operator.

Note that if an interface common to both Petsc and regular matrices is used when accessing the Ap field, it is possible to write functions that operate on both regular and Petsc linear operators.

Public Fields

A: a PetscMat

LinearSolvers.PetscMatLO — Method.

Outer constructor for PetscMatLO

Inputs

pc: a Petsc preconditioner of some kind
mesh
sbp
eqn
opts

LinearSolvers.PetscMatFreeLO — Type.

Petsc matrix-free linear operator.

Public Fields

none

LinearSolvers.PetscMatFreeLO — Method.

Outer constructor for PetscMatFreeLO

Inputs

pc: a Petsc preconditioner of some kind
mesh
sbp
eqn
opts

PetscMatFreePC has an additional function not needed by other linear operator types:

LinearSolvers.setLOCtx — Function.

This function sets the ctx for the underlying Petsc matrix object. The end result is that the mesh, sbp, eqn, and opts that are passed into this function will be passed to the user-defined applyLinearOperator during the next linear solve. See that function for a description of the arguments.

Inputs

lo: the user defined linear operator
mesh
sbp
eqn
opts
ctx_residual
t

Implementing a New Linear Operator

Many of the API functions are automatically generated for matrix-explicit linear operators. Matrix-free linear operators need to define most of the LO interface themselves. A summary of the functions each preconditioner must implement is:

Matrix-explicit

calcLinearOperator

Matrix-free

Example

Linear operators use the same composition structure as preconditioners, see the Implementing a New Preconditioner section before proceeding.

# Assume the ExamplePetscMatLO type has already been defined

type OuterPetscMatLO <: AbstractPetscMatPC
  lo_inner::ExamplePetscMatLO
end


function OuterPetscMatPC(pc::OuterPetscMatLO, mesh::AbstractMesh,
                         sbp::AbstractSBP, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = ExamplePetscMatLO(pc, mesh, sbp, eqn, opts)
  return OuterPetscMatPC(lo_inner)
end

function calcLinearOperator(lo::OuterPetscMatLO, mesh::AbstractMesh,
                            sbp::AbstractSBP, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  # have inner LO do its computation
  calcLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t)

  lo2 = getBaseLO(lo)
  # use set_values1!(lo2.A, ...) to modify the matrix A

  return nothing
end

applyLinearOperator and applyLinearOperatorTranspose are defined automatically for all AbstractPetscMatLO types (using getBaseLO).