Why do we need our wrapper of a pardiso solver?

[samcunliffe]

In solve_*d.py we have a wrapper class for a pardiso solver. Want to understand why we don't simply use an instance of the pardiso solver directly.

Our solve method seems to overlap with the upstream class implementation.

class PySolver:
    def __init__(self, Asp, psolver):  # noqa: N803 | convention Ax = b
        self.Asp = Asp
        self.solver = psolver
        if not pypardiso:
            warnings.warn(
                "Initialising a PySolver, but PyPardiso is not available.",
                stacklevel=2,
            )

    def solve(self, b_inp, x_out):
        self.solver.check_A(self.Asp)
        b = self.solver._check_b(self.Asp, b_inp.array)  # noqa: SLF001, TODO: fix this.
        self.solver.set_phase(33)
        x_out.array[:] = self.solver._call_pardiso(self.Asp, b)[:]  # noqa: SLF001, TODO: fix this.

[janoschpreuss]

This is a hack I introduced to improve performance. During the GMRES iteration the same factorized matrix is used many times. In this case, the original solve method of the pypardiso solver will first perform a check (in every GMRES iteration ...) whether the matrix is already factorized and then set self.set_phase(33) and perform the solver. This check is extremely slow (as it converts things to numpy arrays) and unnecessary, because we know that the matrix is already factorized (our code does this explicitly before starting the GMRes iteration).