def _solve(self, A=None, b=None):
        r"""
        Sends the A and b matrices to the specified solver, and solves for *x*
        given the boundary conditions, and source terms based on the present
        value of *x*.  This method does NOT iterate to solve for non-linear
        source terms or march time steps.

        Parameters
        ----------
        A : sparse matrix
            The coefficient matrix in sparse format. If not specified, then
            it uses  the ``A`` matrix attached to the object.

        b : ND-array
            The RHS matrix in any format.  If not specified, then it uses
            the ``b`` matrix attached to the object.

        Notes
        -----
        The solver used here is specified in the ``settings`` attribute of the
        algorithm.

        """
        # Fetch A and b from self if not given, and throw error if they've not
        # been calculated
        if A is None:
            A = self.A
            if A is None:
                raise Exception('The A matrix has not been built yet')
        if b is None:
            b = self.b
            if b is None:
                raise Exception('The b matrix has not been built yet')
        A = A.tocsr()

        # Default behavior -> use Scipy's default solver (spsolve)
        if self.settings['solver'] == 'pyamg':
            self.settings['solver_family'] = 'pyamg'
        if self.settings['solver'] == 'petsc':
            self.settings['solver_family'] = 'petsc'

        # Set tolerance for iterative solvers
        rtol = self.settings['solver_rtol']
        min_A = np.abs(A.data).min()
        min_b = np.abs(b).min() or 1e100
        atol = min(min_A, min_b) * rtol

        # SciPy
        if self.settings['solver_family'] == 'scipy':
            if importlib.util.find_spec('scikit-umfpack'):
                A.indices = A.indices.astype(np.int64)
                A.indptr = A.indptr.astype(np.int64)
            iterative = ['bicg', 'bicgstab', 'cg', 'cgs', 'gmres', 'lgmres',
                         'minres', 'gcrotmk', 'qmr']
            solver = getattr(sprs.linalg, self.settings['solver_type'])
            if self.settings['solver_type'] in iterative:
                x, exit_code = solver(A=A, b=b, atol=atol, tol=rtol,
                                      maxiter=self.settings['solver_maxiter'])
                if exit_code > 0:
                    raise Exception('SciPy solver did not converge! ' +
                                    'Exit code: ' + str(exit_code))
            else:
                x = solver(A=A, b=b)
            return x

        # PETSc
        if self.settings['solver_family'] == 'petsc':
            # Check if petsc is available
            if importlib.util.find_spec('petsc4py'):
                from openpnm.utils.petsc import PETScSparseLinearSolver as SLS
            else:
                raise Exception('PETSc is not installed.')
            # Define the petsc linear system converting the scipy objects
            ls = SLS(A=A, b=b)
            sets = self.settings
            sets = {k: v for k, v in sets.items() if k.startswith('solver_')}
            sets = {k.split('solver_')[1]: v for k, v in sets.items()}
            ls.settings.update(sets)
            x = SLS.solve(ls)
            del(ls)
            return x

        # PyAMG
        if self.settings['solver_family'] == 'pyamg':
            if importlib.util.find_spec('pyamg'):
                import pyamg
            else:
                raise Exception('pyamg is not installed.')
            ml = pyamg.ruge_stuben_solver(A)
            x = ml.solve(b=b, tol=1e-6)
            return x