Utils

Mat

resolvent4py.utils.matrix.assemble_harmonic_resolvent_generator(A: Mat, freqs: array) → Mat

Assemble \(T = -M + A\), where \(A\) is the output of resolvent4py.utils.io.read_harmonic_balanced_matrix() and \(M\) is a block diagonal matrix with block \(k\) given by \(M_k = i k \omega I\) and \(k\omega\) is the \(k\) th entry of freqs.

Parameters:

• A (PETSc.Mat) – assembled PETSc matrix

• freqs (np.array) – array \(\omega\left(\ldots, -1, 0, 1, \ldots\right)\)

Return type:

PETSc.Mat

resolvent4py.utils.matrix.convert_coo_to_csr(arrays: tuple[array, array, array], sizes: tuple[tuple[int, int], tuple[int, int]]) → tuple[array, array, array]

Convert arrays = [row indices, col indices, values] for COO matrix assembly to [row pointers, col indices, values] for CSR matrix assembly. (Petsc4py currently does not support COO matrix assembly, hence the need to convert.)

Parameters:

• arrays – a list of numpy arrays (e.g., arrays = [rows,cols,vals])

• sizes (tuple[np.array, np.array, np.array]) – see MatSizeSpec

Returns:

csr row pointers, column indices and matrix values for CSR matrix assembly

Return type:

tuple[np.array, np.array, np.array]

resolvent4py.utils.matrix.create_AIJ_identity(comm: Comm, sizes: tuple[tuple[int, int], tuple[int, int]]) → Mat

Create identity matrix of sparse AIJ type

Parameters:

• comm (PETSc.Comm) – MPI Communicator

• sizes (tuple[tuple[int, int], tuple[int, int]]) – see MatSizeSpec

Returns:

identity matrix

Return type:

PETSc.Mat.Type.AIJ

resolvent4py.utils.matrix.create_dense_matrix(comm: Comm, sizes: tuple[tuple[int, int], tuple[int, int]]) → Mat

Create dense matrix

Parameters:

• comm (PETSc.Comm) – PETSc communicator

• sizes – tuple[tuple[int, int], tuple[int, int]]

Return type:

PETSc.Mat.Type.DENSE

resolvent4py.utils.matrix.hermitian_transpose(Mat: Mat, in_place=False, MatHT=None) → Mat

Return the hermitian transpose of the matrix Mat.

Parameters:

• Mat (PETSc.Mat) – PETSc matrix

• in_place (Optional[bool] defaults to False) – in-place transposition if True and out of place otherwise

• MatHT – [optional] matrix with the correct layout to hold the hermitian transpose of Mat

• MatHT – Optional[PETSc.Mat] defaults to None

resolvent4py.utils.matrix.mat_solve_hermitian_transpose(ksp: KSP, X: Mat, Y: Mat | None = None) → Mat

Solve \(A^{-*}X = Y\), where \(X\) is a PETSc matrix of type PETSc.Mat.Type.DENSE

Parameters:

• ksp (PETSc.KSP) – a KPS solver structure

• X (PETSc.Mat.Type.DENSE) – a dense PETSc matrix

• Y (Optional[PETSc.Mat.Type.DENSE] defaults to None) – a dense PETSc matrix

Returns:

matrix to store the result

Return type:

PETSc.Mat.Type.DENSE

Vec

resolvent4py.utils.vector.check_complex_conjugacy(comm: Comm, vec: Vec, nblocks: int) → bool

Verify whether the components \(v_i\) of the vector

\[v = \left(\ldots,v_{-1},v_{0},v_{1},\ldots\right)\]

satisfy \(v_{-i} = \overline{v_{i}}\) for all \(i\).

Parameters:

• vec (PETSc.Vec) – vector \(v\) described above

• nblocks (int) – number of vectors \(v_i\) in \(v\). This must be an odd number.

Returns:

True if the components are complex-conjugates of each other and False otherwise

Return type:

Bool

resolvent4py.utils.vector.enforce_complex_conjugacy(comm: Comm, vec: Vec, nblocks: int) → None

Suppose we have a vector

\[v = \left(\ldots,v_{-1},v_{0},v_{1},\ldots\right)\]

where \(v_i\) are complex vectors. This function enforces \(v_{-i} = \overline{v_{i}}\) for all \(i\) (this implies that \(v_0\) will be purely real).

Parameters:

• vec (PETSc.Vec) – vector \(v\) described above

• nblocks (int) – number of vectors \(v_i\) in \(v\). This must be an odd number.

Return type:

None

resolvent4py.utils.vector.vec_imag(x: Vec, inplace: bool | None = False) → Vec

Returns the imaginary part \(\text{Im}(x)\) of the PETSc Vec.

Parameters:

inplace (Optional[bool], default is False) – in-place if True, else the result is stored in a new PETSc.Vec

Return type:

PETSc.Vec

resolvent4py.utils.vector.vec_real(x: Vec, inplace: bool | None = False) → Vec

Returns the real part \(\text{Re}(x)\) of the PETSc Vec.

Parameters:

inplace (Optional[bool], default is False) – in-place if True, else the result is stored in a new PETSc.Vec

Return type:

PETSc.Vec

BV

resolvent4py.utils.bv.bv_add(alpha: float, X: BV, Y: BV) → None

Compute in-place addition \(X \leftarrow X + \alpha Y\)

Return type:

SLEPc.BV

resolvent4py.utils.bv.bv_conj(X: BV, inplace: bool | None = False) → BV

Returns the complex conjugate \(\overline{X}\) of the BV structure

Parameters:

inplace (Optional[bool], default is False) – in-place if True, else the result is stored in a new SLEPc.BV structure

Return type:

SLEPc.BV

resolvent4py.utils.bv.bv_imag(X: BV, inplace: bool | None = False) → BV

Returns the imaginary part \(\text{Im}(X)\) of the BV structure

Parameters:

inplace (Optional[bool], default is False) – in-place if True, else the result is stored in a new SLEPc.BV structure

Return type:

SLEPc.BV

resolvent4py.utils.bv.bv_real(X: BV, inplace: bool | None = False) → BV

Returns the real part \(\text{Re}(X)\) of the BV structure

Parameters:

inplace (Optional[bool], default is False) – in-place if True, else the result is stored in a new SLEPc.BV structure

Return type:

SLEPc.BV

Ksp

resolvent4py.utils.ksp.check_gmres_bjacobi_solver(A: Mat, ksp: KSP) → None

Check that the solver computed in create_gmres_bjacobi_solver() has succeeded.

Parameters:

• A (PETSc.Mat) – PETSc matrix

• ksp (PETSc.KSP) – PETSc KSP solver

resolvent4py.utils.ksp.check_lu_factorization(A: Mat, ksp: KSP) → None

Check that the LU factorization computed in create_mumps_solver() has succeeded.

Parameters:

• A (PETSc.Mat) – PETSc matrix

• ksp (PETSc.KSP) – PETSc KSP solver

resolvent4py.utils.ksp.create_gmres_bjacobi_solver(A: Mat, nblocks: int, rtol: float | None = 1e-10, atol: float | None = 1e-10, monitor: bool | None = False) → KSP

Create GMRES solver with block-jacobi preconditioner.

Parameters:

• A (PETSc.Mat) – PETSc matrix

• nblocks (int) – number of blocks for the block jacobi preconditioner

• rtol (Optional[float], default is \(10^{-10}\)) – relative tolerance for GMRES

• atol (Optional[float], default is \(10^{-10}\)) – absolute tolerance for GMRES

• monitor (Optional[bool], default is False) – True to monitor convergence and print residual history to terminal. False otherwise

Return ksp:

PETSc KSP solver

Rtype ksp:

PETSc.KSP

resolvent4py.utils.ksp.create_mumps_solver(A: Mat) → KSP

Compute an LU factorization of the matrix A using MUMPS.

Parameters:

A (PETSc.Mat) – PETSc matrix

Return ksp:

PETSc KSP solver

Rtype ksp:

PETSc.KSP

Comms

resolvent4py.utils.comms.compute_local_size(Ng: int) → int

Given the global size Nglob of a vector, compute the local size Nloc that will be owned by each processor in the MPI pool according to the formula

\[\begin{split}N_l = \begin{cases} \left\lfloor \frac{N_g}{S} \right\rfloor + 1 & \text{if} \, N_g \text{mod} S > r \\ \left\lfloor \frac{N_g}{S} \right\rfloor & \text{otherwise} \end{cases}\end{split}\]

where \(N_g\) is the global size, \(N_l\) is the local size, \(S\) is the size of the MPI pool (i.e., the total number of processors) and \(r\) is the rank of the current processor.

Parameters:

Ng (int) – global size

Returns:

local size

Return type:

int

resolvent4py.utils.comms.distributed_to_sequential_matrix(Mat_dist: Mat) → Mat

Allgather dense distributed PETSc matrix into a dense sequential PETSc matrix.

Parameters:

Mat_dist (PETSc.Mat.Type.DENSE) – a distributed dense matrix

Return type:

PETSc.Mat.Type.DENSE

resolvent4py.utils.comms.distributed_to_sequential_vector(vec_dist: Vec) → Vec

Allgather distribued PETSc vector into PETSc sequential vector.

Parameters:

vec_dist (PETSc.Vec.Type.STANDARD) – a distributed vector

Return type:

PETSc.Vec.Type.SEQ

resolvent4py.utils.comms.scatter_array_from_root_to_all(array: array, locsize: int | None = None) → array

Scatter numpy array from root to all other processors

Parameters:

• array (numpy.array) – numpy array to be scattered

• locsize (Optional[int] default to None) – local size owned by each processor. If None, this is computed using the same logic as in the compute_local_size() routine.

Return type:

numpy.array

resolvent4py.utils.comms.sequential_to_distributed_matrix(Mat_seq: Mat, Mat_dist: Mat) → Mat

Scatter a sequential dense PETSc matrix (type PETSc.Mat.Type.DENSE) to a distributed dense PETSc matrix (type PETSc.Mat.Type.DENSE)

Parameters:

• Mat_seq (PETSc.Mat.Type.DENSE) – a sequential dense matrix

• Mat_dist (PETSc.Mat.Type.DENSE) – a distributed dense matrix

Returns:

a distributed dense matrix

Return type:

PETSc.Mat.Type.DENSE

resolvent4py.utils.comms.sequential_to_distributed_vector(vec_seq: Vec, vec_dist: Vec) → Vec

Scatter a sequential PETSc vector (type PETSc.Vec.Type.SEQ) to a distributed PETSc vector (type PETSc.Vec.Type.STANDARD)

Parameters:

• vec_seq (PETSc.Vec.Type.SEQ) – a sequential vector

• vec_dist (PETSc.Vec.Type.STANDARD) – a distributed vector

Returns:

a distributed vector

Return type:

PETSc.Vec.Type.STANDARD

I/O

resolvent4py.utils.io.read_bv(filename: str, sizes: Tuple[Tuple[int, int], int]) → BV

Read dense matrix from file and store as a SLEPc BV

Parameters:

• filename (str) – name of the file that holds the matrix

• sizes (Tuple[Tuple[int, int], int]) – ((local rows, global rows), global columns)

Return type:

SLEPc.BV

resolvent4py.utils.io.read_coo_matrix(filenames: Tuple[str, str, str], sizes: Tuple[Tuple[int, int], Tuple[int, int]]) → Mat

Read COO matrix from file

Parameters:

• filenames (Tuple[str, str, str]) – names of the files that hold the rows, columns and values of the sparse matrix (e.g., (rows, cols, values))

• sizes (Tuple[Tuple[int, int], Tuple[int, int]]) – ((local rows, global rows), (local cols, global cols))

Return type:

PETSc.Mat

resolvent4py.utils.io.read_dense_matrix(filename: str, sizes: Tuple[Tuple[int, int], Tuple[int, int]]) → Mat

Read dense PETSc matrix from file.

Parameters:

• filename (str) – name of the file that holds the matrix

• sizes (Tuple[Tuple[int, int], Tuple[int, int]]) – ((local rows, global rows), (local cols, global cols))

Return type:

PETSc.Mat

resolvent4py.utils.io.read_harmonic_balanced_bv(filenames_lst: List[str], real_bflow: bool, block_sizes: Tuple[Tuple[int, int], int], full_sizes: Tuple[Tuple[int, int], int]) → BV

Given \(\{\ldots, A_{-1}, A_{0}, A_{1},\ldots\}\), where \(A_j\) is the \(j\) th Fourier coefficient of a time-periodic matrix \(A(t)\), assemble the harmonic-balanced matrix

\[\begin{split}M = \begin{bmatrix} \ddots & \ddots & \ddots & \ddots \\ \ddots & A_{0} & A_{-1} & A_{-2} & \ddots \\ \ddots & A_{1} & A_{0} & A_{-1} & A_{-2} & \ddots \\ \ddots & A_{2} & A_{1} & A_{0} & A_{-1} & A_{-2} &\ddots \\ & \ddots & A_{2} & A_{1} & A_{0} & A_{-1} &\ddots \\ & & \ddots & A_{2} & A_{1} & A_{0} & \ddots\\ & & & \ddots & \ddots & \ddots & \ddots \\ \end{bmatrix}\end{split}\]

Parameters:

• comm (PETSc.Comm) – MPI communicator (PETSc.COMM_WORLD)

• filenames_lst (List[str]) – list of names of files where the Fourier modes \(A_j\) are stored as SLEPc BVs

• real_bflow (bool) – set to True if filenames_lst contains the names of the COO arrays of the positive Fourier coefficients of \(A(t)\) (i.e., \(\{A_0, A_1, A_2, \ldots\}\)). The negative frequencies are assumed the complex-conjugates of the positive ones. Set to False otherwise (i.e., \(\{\ldots, A_{-2}, A_{-1}, A_0, A_1, A_2, \ldots\}\)).

• block_sizes (Tuple[Tuple[int, int], int]) – sizes ((local rows, global rows), global_cols) of \(A_j\)

• full_sizes (Tuple[Tuple[int, int], int]) – sizes ((local rows, global rows), global_cols) of \(M\)

Returns:

SLEPc BV \(M\)

Return type:

SLEPc.BV

resolvent4py.utils.io.read_harmonic_balanced_matrix(filenames_lst: List[Tuple[str, str, str]], real_bflow: bool, block_sizes: Tuple[Tuple[int, int], Tuple[int, int]], full_sizes: Tuple[Tuple[int, int], Tuple[int, int]]) → Mat

Given \(\{\ldots, A_{-1}, A_{0}, A_{1},\ldots\}\), where \(A_j\) is the \(j\) th Fourier coefficient of a time-periodic matrix \(A(t)\), assemble the harmonic-balanced matrix

\[\begin{split}M = \begin{bmatrix} \ddots & \ddots & \ddots & \ddots \\ \ddots & A_{0} & A_{-1} & A_{-2} & \ddots \\ \ddots & A_{1} & A_{0} & A_{-1} & A_{-2} & \ddots \\ \ddots & A_{2} & A_{1} & A_{0} & A_{-1} & A_{-2} &\ddots \\ & \ddots & A_{2} & A_{1} & A_{0} & A_{-1} &\ddots \\ & & \ddots & A_{2} & A_{1} & A_{0} & \ddots\\ & & & \ddots & \ddots & \ddots & \ddots \\ \end{bmatrix}\end{split}\]

Parameters:

• filenames_lst (List[Tuple[str, str, str]]) – list of tuples, with each tuple of the form (rows_j.dat, cols_j.dat, vals_j.dat) containing the COO arrays of the matrix \(A_j\)

• real_bflow (bool) – set to True if filenames_lst contains the names of the COO arrays of the positive Fourier coefficients of \(A(t)\) (i.e., \(\{A_0, A_1, A_2, \ldots\}\)). The negative frequencies are assumed the complex-conjugates of the positive ones. Set to False otherwise (i.e., \(\{\ldots, A_{-2}, A_{-1}, A_0, A_1, A_2, \ldots\}\)).

• block_sizes (Tuple[Tuple[int, int], Tuple[int, int]]) – sizes ((local rows, global rows), (local cols, global cols)) of \(A_j\)

• full_sizes (Tuple[Tuple[int, int], Tuple[int, int]]) – sizes ((local rows, global rows), (local cols, global cols)) of \(M\)

Returns:

Sparse PETSc matrix \(M\)

Return type:

PETSc.Mat

resolvent4py.utils.io.read_harmonic_balanced_vector(filenames_lst: List[str], real_bflow: bool, block_sizes: Tuple[int, int], full_sizes: Tuple[int, int]) → Vec

Given \(\{\ldots, v_{-1}, v_{0}, v_{1},\ldots\}\), where \(v_j\) is the \(j\) th Fourier coefficient of a time-periodic vector \(v(t)\), assemble

\[\begin{split}\hat{v} = \begin{bmatrix} \vdots \\ v_{-1} \\ v_0 \\ v_{1} \\ \vdots \end{bmatrix}.\end{split}\]

Parameters:

• filenames_lst (List[str]) – list of names of files where the Fourier modes \(v_j\) are stored as PETSc Vec

• real_bflow (bool) – set to True if filenames_lst contains the names of the COO arrays of the positive Fourier coefficients of \(v(t)\) (i.e., \(\{v_0, v_1, v_2, \ldots\}\)). The negative frequencies are assumed the complex-conjugates of the positive ones. Set to False otherwise (i.e., \(\{\ldots, v_{-2}, v_{-1}, v_0, v_1, v_2, \ldots\}\)).

• block_sizes (Tuple[int, int]) – sizes (local rows, global rows) of \(v_j\)

• full_sizes (Tuple[int, int]) – sizes (local rows, global rows) of \(\hat{v}\)

Returns:

vector \(\hat{v}\)

Return type:

PETSc.Vec

resolvent4py.utils.io.read_vector(filename: str, sizes: Tuple[int, int] | None = None) → Vec

Read PETSc vector from file

Parameters:

• filename (str) – name of the file that holds the vector

• sizes (Optional[Tuple[int, int]]) – (local size, global size)

Return type:

PETSc.Vec

resolvent4py.utils.io.write_to_file(filename: str, object: Mat | Vec | BV) → None

Write PETSc/SLEPc object (PETSc.Mat, PETSc.Vec or SLEPc.BV) to file.

Parameters:

• filename (str) – name of the file to store the object

• object (Union[PETSc.Mat, PETSc.Vec, SLEPc.BV]) – any PETSc matrix or vector, or SLEPc BV

Random

resolvent4py.utils.random.generate_random_petsc_sparse_matrix(sizes: Tuple[Tuple[int, int], Tuple[int, int]], nnz: int, complex: bool | None = False) → Mat

Parameters:

• sizes (Tuple[Tuple[int, int], Tuple[int, int]]) – ((local rows, global rows), (local cols, global cols))

• nnz (int) – number of non-zero entries

• complex (Optional[bool], default is False) – True if you want a complex-valued matrix

Returns:

a sparse PETSc matrix

Return type:

PETSc.Mat

resolvent4py.utils.random.generate_random_petsc_vector(sizes: Tuple[int, int], complex: bool | None = False) → Vec

Parameters:

• sizes (Tuple[int, int]) – vector sizes

• complex (Optional[bool], default is False) – True if you want a complex-valued vector

Return type:

PETSc.Vec

Miscellaneous

resolvent4py.utils.miscellaneous.get_memory_usage(comm: Comm) → float

Compute the used memory (in Mb) across the MPI pool

Return type:

float

resolvent4py.utils.miscellaneous.get_mpi_type(dtype: dtype) → Datatype

Get the corresponding MPI type for a given numpy data type.

Parameters:

dtype (np.dtypes) – (e.g., np.dtype(PETSc.IntType))

Return type:

MPI.Datatype

resolvent4py.utils.miscellaneous.petscprint(comm: Comm, arg: any) → None

Print to terminal

Parameters:

• comm (PETSc.Comm) – MPI communicator (PETSc.COMM_WORLD or PETSc.COMM_SELF)

• arg (any) – argument to be fed into print()