Linear Operators

Linear Operator

class resolvent4py.linear_operators.linear_operator.LinearOperator(comm: Comm, name: str, dimensions: tuple[tuple[int, int], tuple[int, int]], nblocks: int | None = None)[source]

Bases: object

Abstract base class for linear operators \(L\).

Parameters:

comm (PETSc.Comm) – MPI communicator PETSc.COMM_WORLD
name (str) – name of the linear operator
dimensions (tuple[tuple[int, int], tuple[int, int]]) – local and global sizes of the range and domain of the linear operator, [[local rows, global rows], [local columns, global columns]]
nblocks (Optional[Unions[int, None]], default is None) – number of blocks (if the linear operator has block structure)

abstract apply(x: Vec, y: Vec | None = None) → Vec[source]

Compute \(y = Lx\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose(x: Vec, y: Vec | None = None) → Vec[source]

Compute \(y = L^*x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose_mat(X: BV, Y: BV | None = None) → BV[source]

Compute \(Y = L^*X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

abstract apply_mat(X: BV, Y: BV | None = None) → BV[source]

Compute \(Y = LX\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

check_if_complex_conjugate_structure() → bool[source]

Given a vector

\[x = \left(\ldots,x_{-1},x_{0},x_{1},\ldots\right)\]

with vector-valued entries that satisfy \(x_{-i} = \overline{x_i}\), check if the vector \(Lx\) satisfies \((Lx)_{-i}=\overline{(Lx)_{i}}\). (Here, the overline denote complex conjugation.)

Returns:: True if the linear operator has complex-conjugate structure, False otherwise.
Return type:: bool

check_if_real_valued() → bool[source]

Returns:: True if the linear operator is real-valued, False otherwise
Return type:: bool

create_left_bv(ncols: int) → BV[source]

Parameters:

ncols – number of columns in the BV
type – int

Returns:

a SLEPc BV where \(LX\) can be stored into

Return type:

SLEPc.BV

create_left_vector() → Vec[source]

Returns:: a PETSc vector where \(Lx\) can be stored into
Return type:: PETSc.Vec

create_right_bv(ncols: int) → BV[source]

Parameters:: ncols (int) – number of columns in the BV
Returns:: a SLEPc BV that \(L\) can be multiplied against
Return type:: SLEPc.BV

create_right_vector() → Vec[source]

Returns:: a PETSc vector that \(L\) can be multiplied against
Return type:: PETSc.Vec

abstract destroy() → None[source]: Destroy the PETSc and SLEPc objects associated with \(L\)

get_block_cc_flag() → bool | None[source]

Returns True if the operator has complex-conjugate block structure (like in the harmonic-balancing operator). False if the operator does not have complex-conjugate block structure, None if .get_nblocks returns None.

Return type:: Union[bool, None]

get_comm() → Comm[source]

The MPI communicator

Return type:: PETSc.Comm

get_dimensions() → tuple[tuple[int, int], tuple[int, int]][source]

The local and global dimensions of the linear operator

Return type:: tuple[tuple[int, int], tuple[int, int]]

get_name() → str[source]

The name of the linear operator

Return type:: str

get_nblocks() → int | None[source]

The number of blocks of the linear operator

Return type:: Union[int, None]

get_real_flag() → bool[source]

True if the operator is real-valued, False otherwise.

Return type:: bool

solve(x: Vec, y: Vec | None = None) → Vec[source]

Compute \(y = L^{-1}x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

solve_hermitian_transpose(x: Vec, y: Vec | None = None) → Vec[source]

Compute \(y = L^{-*}x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

solve_hermitian_transpose_mat(X: BV, Y: BV | None = None) → BV[source]

Compute \(Y = L^{-*}X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

solve_mat(X: BV, Y: BV | None = None) → BV[source]

Compute \(Y = L^{-1}X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

Matrix Linear Operator

class resolvent4py.linear_operators.matrix.MatrixLinearOperator(A: Mat, ksp: KSP | None = None, nblocks: int | None = None)[source]

Bases: LinearOperator

Class for a linear operator \(L = A\), where \(A\) is a matrix. In general, A can be any matrix (rectangular or square, invertible or non-invertible). For an invertible \(A\), the user may also provide a PETSc KSP object ksp to act with \(A^{-1}\) on vectors and matrices.

Parameters:

A (PETSc.Mat) – a PETSc matrix
ksp (Optional[Union[PETSc.KSP, None]], defaults to None) – a PETSc KSP object to enable the solve() and solve_hermitian_transpose() methods
nblocks (Optional[Union[int, None]], default is None) – number of blocks (if the linear operator has block structure). This must be an odd number.

apply(x, y=None)[source]

Compute \(y = Lx\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose(x, y=None)[source]

Compute \(y = L^*x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose_mat(X, Y=None)[source]

Compute \(Y = L^*X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

apply_mat(X, Y=None)[source]

Compute \(Y = LX\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

destroy()[source]: Destroy the PETSc and SLEPc objects associated with \(L\)

solve(x, y=None)[source]

Compute \(y = L^{-1}x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

solve_hermitian_transpose(x, y=None)[source]

Compute \(y = L^{-*}x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

solve_hermitian_transpose_mat(X, Y=None)[source]

Compute \(Y = L^{-*}X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

solve_mat(X, Y=None)[source]

Compute \(Y = L^{-1}X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

Low-Rank Linear Operator

class resolvent4py.linear_operators.low_rank.LowRankLinearOperator(U: BV, Sigma: ndarray, V: BV, nblocks: int | None = None)[source]

Bases: LinearOperator

Class for a linear operator of the form

\[L = U \Sigma V^*,\]

where \(U\), \(\Sigma\) and \(V\) are matrices of conformal sizes (and \(\Sigma\) is not necessarily diagonal).

Parameters:

U (SLEPc.BV) – a tall and skinny matrix
Sigma (numpy.ndarray) – 2D numpy array
V (SLEPc.BV) – a tall and skinny matrix
nblocks (Optional[Union[int, None]], default is None) – number of blocks (if the operator has block structure)

apply(x, y=None)[source]

Compute \(y = Lx\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose(x, y=None)[source]

Compute \(y = L^*x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose_mat(X, Y=None)[source]

Compute \(Y = L^*X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

apply_mat(X, Y=None)[source]

Compute \(Y = LX\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

destroy()[source]: Destroy the PETSc and SLEPc objects associated with \(L\)

Low-Rank-Updated Linear Operator

class resolvent4py.linear_operators.low_rank_updated.LowRankUpdatedLinearOperator(A: LinearOperator, B: BV, K: ndarray, C: BV, woodbury_factors: Tuple[BV, ndarray, BV] | None = None, nblocks: int | None = None)[source]

Bases: LinearOperator

Class for a linear operator of the form

\[L = A + B K C^*\]

where \(A\) is an instance of the LinearOperator class, and \(B\), \(K\) and \(C\) are low-rank (dense) matrices of conformal sizes. If A.solve() is enabled, then the solve() method in this class is implemented using the Woodbury matrix identity

\[L^{-1} = A^{-1} - X D Y^*,\]

where \(X\), \(D\) and \(Y\) the Woodbury factors defined as

\[\textcolor{black}{X} = A^{-1}B,\quad \textcolor{black}{D} = K\left(I + C^* A^{-1}B K\right)^{-1},\quad \textcolor{black}{Y} = A^{-*}C,\]

Parameters:

A – instance of the LinearOperator class
B (SLEPc.BV) – tall and skinny matrix
K (numpy.ndarray) – dense matrix
C (SLEPc.BV) – tall and skinny matrix
woodbury_factors (Optional[Union[Tuple[SLEPc.BV, numpy.ndarray, SLEPc.BV], None]], default is None) – tuple \((X, D, Y)\) of Woodbury factors. If A.solve() is enabled and the argument woodbury_factors is None, the factors \(X\), \(D\) and \(Y\) are computed at initialization
nblocks (Optional[Union[int, None]], default is None) – number of blocks (if the operator has block structure)

apply(x, y=None)[source]

Compute \(y = Lx\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose(x, y=None)[source]

Compute \(y = L^*x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose_mat(X, Y=None, Z=None)[source]

Compute \(Y = L^*X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

apply_mat(X, Y=None, Z=None)[source]

Compute \(Y = LX\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

compute_woodbury_operator(nblocks: int | None) → LowRankLinearOperator[source]

Parameters:: nblocks (Unions[int, None]) – number of blocks (if the operator has block structure)
Returns:: a LowRankLinearOperator constructed from the Woodbury factors X, D and Y
Return type:: LowRankUpdatedLinearOperator

create_intermediate_bv(m: int) → BV[source]

Create matrix to store \(X D Y^* Z\), for any matrix Z with \(m\) columns.

Attention

It is the user’s responsibility to destroy this object when no longer needed.

Parameters:: m (int) – number of columns of \(Z\)
Return type:: SLEPc.BV

create_intermediate_bv_hermitian_transpose(m: int) → BV[source]

Create matrix to store \(Y D X^* Z\), for any matrix Z with \(m\) columns.

Attention

It is the user’s responsibility to destroy this object when no longer needed.

Parameters:: m (int) – number of columns of \(Z\)
Return type:: SLEPc.BV

destroy()[source]: Destroy the PETSc and SLEPc objects associated with \(L\)

solve(x, y=None)[source]

Compute \(y = L^{-1}x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

solve_hermitian_transpose(x, y=None)[source]

Compute \(y = L^{-*}x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

solve_hermitian_transpose_mat(X, Y=None, Z=None)[source]

Compute \(Y = L^{-*}X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

solve_mat(X, Y=None, Z=None)[source]

Compute \(Y = L^{-1}X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

Product Linear Operator

class resolvent4py.linear_operators.product.ProductLinearOperator(linops: List[LinearOperator], linops_actions: List[Callable[[Vec, Vec | None], Vec]], nblocks: int | None = None)[source]

Bases: LinearOperator

Linear operator of the form

\[L = L_{r}L_{r-1}\ldots L_{2}L_{1}\]

where \(L_i\) are instances of the LinearOperator class. For example, by setting linops = (L2, L1) and linops_actions = (L2.solve, L1.apply_hermitian_transpose) we define a linear operator

\[L = L_2^{-1}L_1^*.\]

Parameters:

linops (List[LinearOperator]) – list of linear operators (see example above).
linops_actions (List[Callable[[PETSc.Vec, PETSc.Vec], PETSc.Vec]]) – list of actions (see example above). Must be one of LinearOperator.apply(), LinearOperator.solve(), LinearOperator.apply_hermitian_transpose() or LinearOperator.solve_hermitian_transpose()
nblocks (Optional[Union[int, None]], default is None) – number of blocks (if the operator has block structure)

apply(x, y=None)[source]

Compute \(y = Lx\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose(x, y=None)[source]

Compute \(y = L^*x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose_mat(X, Y=None, Z=None)[source]

Compute \(Y = L^*X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

apply_mat(X, Y=None, Z=None)[source]

Compute \(Y = LX\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

create_intermediate_bvs(m: int) → List[BV][source]

Create list of matrices to store \(L_i Z\), for any matrix Z with \(m\) columns.

Attention

It is the user’s responsibility to destroy this object when no longer needed.

Parameters:: m (int) – number of columns of \(Z\)
Return type:: List[SLEPc.BV]

create_intermediate_bvs_hermitian_transpose(m: int) → List[BV][source]

Create list of matrices to store \(L_i^* Z\), for any matrix Z with \(m\) columns.

Attention

It is the user’s responsibility to destroy this object when no longer needed.

Parameters:: m (int) – number of columns of \(Z\)
Return type:: List[SLEPc.BV]

destroy()[source]: Destroy the PETSc and SLEPc objects associated with \(L\)

Projection Linear Operator

class resolvent4py.linear_operators.projection.ProjectionLinearOperator(Phi: BV, Psi: BV, complement: bool | None = False, nblocks: int | None = None)[source]

Bases: LinearOperator

Class for a linear operator of the form

\[\begin{split}L = \begin{cases} \Phi \left(\Psi^*\Phi\right)^{-1} \Psi^{*} & \text{if } \texttt{complement=False} \\ I - \Phi \left(\Psi^*\Phi\right)^{-1} \Psi^{*} & \text{otherwise} \end{cases}\end{split}\]

where \(\Phi\) and \(\Psi\) are tall and skinny matrices of size \(N\times r\) stored as SLEPc BVs. (It is easy to check that \(L\) is a projection since \(L^2 = L\).)

Parameters:

Phi (SLEPc.BV) – tall and skinny matrix
Psi (SLEPc.BV) – tall and skinny matrix
complement (Optional[bool], default if False) – see definition of \(L\)
nblocks (Optional[Union[int, None]], default is None) – number of blocks (if the operator has block structure)

apply(x, y=None)[source]

Compute \(y = Lx\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose(x, y=None)[source]

Compute \(y = L^*x\).

Parameters:

x (PETSc.Vec) – a PETSc vector
y (Optional[PETSc.Vec], default is None) – a PETSc vector to store the result

Return type:

PETSc.Vec

apply_hermitian_transpose_mat(X, Y=None)[source]

Compute \(Y = L^*X\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

apply_mat(X, Y=None)[source]

Compute \(Y = LX\)

Parameters:

X (SLEPc.BV) – a SLEPc BV
Y (Optional[SLEPc.BV], default is None) – a SLEPc BV to store the result

Return type:

SLEPc.BV

destroy()[source]: Destroy the PETSc and SLEPc objects associated with \(L\)