PetscMatrixBase & operator= (PetscMatrixBase &&) = delete;
  virtual ~PetscMatrixBase ();

  virtual SolverPackage solver_package() override
  {
    return PETSC_SOLVERS;
  }

  /**

template <typename T>
void
PetscMFFDMatrix<T>::assign(Mat m, bool set_context)
{
  if (this->_mat != m)
    this->clear();

  this->_mat = m;
  this->_is_initialized = true;
  this->_destroy_mat_on_exit = false;
  if (set_context)
    this->set_context();
}

template <typename T>
PetscMFFDMatrix<T> &

   * the preconditioning matrix -- which may be the same object (the common case) or genuinely
   * distinct (e.g. a matrix-free operator paired with an assembled preconditioning matrix).
   */
  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 */)
  {
    libmesh_not_implemented();
  }

  /**

template <typename T>
void
PetscMatrixShellMatrix<T>::init(const numeric_index_type m,
                                const numeric_index_type n,
                                const numeric_index_type m_l,
                                const numeric_index_type n_l,

                                const numeric_index_type,
                                const numeric_index_type blocksize)
{
  init_shell_mat(*this, m, n, m_l, n_l, blocksize);
  this->set_context();
}

template <typename T>
void

template <typename T>
void
PetscMatrixShellMatrix<T>::zero()
{
  // A shell matrix generally computes its action and stores no entries, so there is nothing to
  // clear. This is reachable through System::init_matrices(), which zeroes every matrix it

  // init_matrices() call with something like an is_shell() attribute, as we consider it relatively
  // harmless to allow a user to "zero" a shell compared to attempting to add/set something
  // nontrivial in the shell
}

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

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_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add(const numeric_index_type, const numeric_index_type, const T)
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

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_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
                                      const std::vector<numeric_index_type> &)
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add(const T, const SparseMatrix<T> &)
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
T
PetscMatrixShellMatrix<T>::operator()(const numeric_index_type, const numeric_index_type) const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
Real
PetscMatrixShellMatrix<T>::l1_norm() const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
Real
PetscMatrixShellMatrix<T>::linfty_norm() const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::print_personal(std::ostream &) const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_diagonal(NumericVector<T> &) const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_transpose(SparseMatrix<T> &) const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_row(numeric_index_type,
                                   std::vector<numeric_index_type> &,
                                   std::vector<T> &) const
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

template <typename T>
SparseMatrix<T> &
PetscMatrixShellMatrix<T>::operator=(const SparseMatrix<T> &)
{
  libmesh_error_msg("Method not appropriate for arbitrary shell matrices");
}

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") ?

  else
    {
      if (this->prefix_with_name())
        nonlinear_solver->init(this->prefix().c_str());
      else
        nonlinear_solver->init();

      // 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