template <typename T>
void
PetscMFFDMatrix<T>::assign(Mat m, bool set_context)
{
  this->_mat = m;
  this->_is_initialized = true;
  this->_destroy_mat_on_exit = false;
  if (set_context)
    this->set_context();
}

template <typename T>
void

   * default implementation ignores \p jac_in and simply forwards to the single-matrix solve(),
   * for solver backends that do not distinguish between the two.
   */
  virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> & /* jac_in */,
                                               SparseMatrix<T> & pre_in,
                                               NumericVector<T> & x_in,
                                               NumericVector<T> & r_in,
                                               const double tol,
                                               const unsigned int m_its)
  {
    return this->solve(pre_in, x_in, r_in, tol, m_its);
  }

  /**

    {
      ray_target = inside - Point(1);
      intersection_distances =
        this->find_ray_intersections(inside, ray_target);
    }

  // I'd make this an assert, but I'm not 100% confident we can't

template <typename T>
void
PetscMatrixShellMatrix<T>::zero()
{
  libmesh_error();
}

template <typename T>
std::unique_ptr<SparseMatrix<T>>
PetscMatrixShellMatrix<T>::zero_clone() const
{
  libmesh_error();
}

template <typename T>
std::unique_ptr<SparseMatrix<T>>
PetscMatrixShellMatrix<T>::clone() const
{
  libmesh_not_implemented();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::set(const numeric_index_type, const numeric_index_type, const T)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add(const numeric_index_type, const numeric_index_type, const T)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
                                      const std::vector<numeric_index_type> &,
                                      const std::vector<numeric_index_type> &)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
                                      const std::vector<numeric_index_type> &)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add(const T, const SparseMatrix<T> &)
{
  libmesh_error();
}

template <typename T>
T
PetscMatrixShellMatrix<T>::operator()(const numeric_index_type, const numeric_index_type) const
{
  libmesh_error();
}

template <typename T>
Real
PetscMatrixShellMatrix<T>::l1_norm() const
{
  libmesh_error();
}

template <typename T>
Real
PetscMatrixShellMatrix<T>::linfty_norm() const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::print_personal(std::ostream &) const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_diagonal(NumericVector<T> &) const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_transpose(SparseMatrix<T> &) const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_row(numeric_index_type,
                                   std::vector<numeric_index_type> &,
                                   std::vector<T> &) const
{
  libmesh_error();
}

template class LIBMESH_EXPORT PetscMatrixShellMatrix<Number>;

        Jac = &mffd_jac;
        // mffd_jac is function-local, so don't attach a context to jac here -- it would
        // dangle once mffd_jac is destroyed at the end of this call.
        mffd_jac.assign(jac, /*set_context=*/false);
      }

    // We already computed the Jacobian during the residual evaluation

  _default_monitor(true),
  _snesmf_reuse_base(true),
  _computing_base_vector(true),
  _setup_reuse(false)
{
}

                                const double        tol,     // Stopping tolerance
                                const unsigned int  m_its)
{
  return this->solve(pre_in, pre_in, x_in, r_in, tol, m_its);
}

template <typename T>
std::pair<unsigned int, Real>
PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  jac_in,  // Jacobian operator matrix (Amat)
                                SparseMatrix<T> &  pre_in,  // Preconditioning matrix (Pmat)
                                NumericVector<T> & x_in,    // Solution vector
                                NumericVector<T> & r_in,    // Residual vector

  this->init ();

  // Make sure the data passed in are really of Petsc types
  PetscMatrixBase<T> * jac = cast_ptr<PetscMatrixBase<T> *>(&jac_in);
  PetscMatrixBase<T> * pre = cast_ptr<PetscMatrixBase<T> *>(&pre_in);
  PetscVector<T> * x   = cast_ptr<PetscVector<T> *>(&x_in);
  PetscVector<T> * r   = cast_ptr<PetscVector<T> *>(&r_in);

  // Only set the jacobian function if we've been provided with something to call.
  // This allows a user to set their own jacobian function if they want to
  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
    LibmeshPetscCall(SNESSetJacobian (_snes, jac->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));

  // Have the Krylov subspace method use our good initial guess rather than 0
  KSP ksp;

  nonlinear_solver          (NonlinearSolver<Number>::build(*this)),
  diff_solver               (),
  _n_nonlinear_iterations   (0),
  _final_nonlinear_residual (1.e20),
  _operator_matrix          (nullptr)
{
  // Set default parameters
  // These were chosen to match the Petsc defaults

    es.parameters.get<unsigned int>("nonlinear solver maximum function evaluations");

  const double abs_resid_tol = parameters.have_parameter<Real>("nonlinear solver absolute residual tolerance") ?
    double(parameters.get<Real>("nonlinear solver absolute residual tolerance")) :
    double(es.parameters.get<Real>("nonlinear solver absolute residual tolerance"));

  const double rel_resid_tol = parameters.have_parameter<Real>("nonlinear solver relative residual tolerance") ?

      // If a distinct Jacobian operator matrix has been registered (see
      // set_operator_matrix()), use it as Amat while *matrix remains the preconditioning matrix
      // (Pmat); otherwise use the ordinary single-matrix solve.
      if (_operator_matrix)
        std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
          nonlinear_solver->solve (*_operator_matrix, *matrix, *solution, *rhs,
                                   nonlinear_solver->relative_residual_tolerance,
                                   nonlinear_solver->max_linear_iterations);
      else
        std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
          nonlinear_solver->solve (*matrix, *solution, *rhs,
                                   nonlinear_solver->relative_residual_tolerance,
                                   nonlinear_solver->max_linear_iterations);
    }

  // Update the system after the solve