pymor.solvers

Additional pyMOR solvers

Imports

Testing

API

diagonal_blocks

def diagonal_blocks(
    mat, block_size
):

source

BlockSpILUPreconditioner

def BlockSpILUPreconditioner(
    matrix:NoneType=None, block_size:NoneType=None, spilu_drop_tol:float=0.0001, spilu_fill_factor:int=20,
    spilu_drop_rule:NoneType=None, spilu_permc_spec:str='COLAMD', n_workers:NoneType=None,
    single_precision:bool=False
):

Equation solver.

Solves operator equations of the form

.. math:: A(U; ) = V

The operator :math:A can be linear or non-linear. When :math:A is linear, a solver can also be used to solve the adjoint equation

.. math:: A^H(V; ) = U

for :math:U.

When :attr:~pymor.solvers.interface.Solver.least_squares is True, the equations are solved in a least-squares sense:

.. math:: {U} |A(U; ) - V|^2 {V} |A^H(V; ) - U|^2

Solvers will typically only work for certain classes of |Operators|. In most cases, solvers are invoked by the :meth:~pymor.operators.interface.Operator.apply_inverse and :meth:~pymor.operators.interface.Operator.apply_inverse_adjoint methods of |Operators|. If an |Operator| has no associated solver, :class:~pymor.solvers.default.DefaultSolver is used.

# from sksparse.umfpack import umf_factor, UMFFactor

# #|export
# from fastcore.parallel import parallel
# from functools import partial
# from sksparse.umfpack import umf_factor, UMFFactor

# class BlockUMFPreconditioner(Solver):
#     def __init__(self, matrix=None, block_size=None):
#         self.matrix = matrix
#         self.block_size = block_size
#         self._ilu = None
#         self._umf = None

#     def factorize(self, mu):
#         if self.matrix is None or self.block_size is None:
#             raise ValueError("Matrix and block size must be specified for the preconditioner.")
#         matrix = to_matrix(self.matrix.assemble(mu)).tocsc()
#         blocks = diagonal_blocks(matrix, self.block_size)
#         self._n_blocks = len(blocks)
#         self._blocks = blocks
#         if self._umf is None: self._umf = UMFFactor(blocks[0])
#         if self._ilu is None: self._ilu = [None] * len(blocks)
#         new_blocks = [i for i in range(len(blocks)) if self._ilu[i] is None]
#         tick = clock()
#         # new_ilus = parallel(
#         #     lambda A: self._umf.copy().factorize(A=A),
#         #     [blocks[i] for i in new_blocks],
#         #     threadpool=True,
#         #     n_workers=4,
#         # )
#         umf = self._umf.copy()
#         umf.factorize(blocks[new_blocks[0]])
#         umf = self._umf.copy()
#         umf.factorize(blocks[new_blocks[1]])
#         new_ilus = [self._umf.copy().factorize(A=blocks[i]) for i in new_blocks]
#         logger.info(f'BlockUMFPreconditioner: Factoring {len(new_blocks)} new blocks out of {len(blocks)} took {clock()-tick:.2f} seconds')
#         for idx, ilu in zip(new_blocks, new_ilus):
#             self._ilu[idx] = ilu

#     def operator(self, mu):
#         self.factorize(mu)
#         ilu = self._ilu
#         self._n_solves = 0

#         def ilu_block_solve(v):
#             v = v.reshape(self._n_blocks, self.block_size)
#             result = np.empty_like(v).astype(self._blocks[0].dtype)
#             for i in range(self._n_blocks):
#                 result[i] = ilu[i].solve(v[i])
#             self._n_solves += 1
#             return result.reshape(-1)
        
#         logger.info(f'BlockUMFPreconditioner: {self._n_blocks} preconditioner blocks of size {self.block_size}, type {self._ilu[0].L.dtype}, and density {self._blocks[0].nnz / (self.block_size**2):.2e}')
#         return LinearOperator((self._n_blocks*self.block_size, self._n_blocks*self.block_size), ilu_block_solve)