diff --git a/Src/Base/AMReX_TableData.H b/Src/Base/AMReX_TableData.H
index 842225e53f4..9d0db5bd8d1 100644
--- a/Src/Base/AMReX_TableData.H
+++ b/Src/Base/AMReX_TableData.H
@@ -22,8 +22,7 @@ struct Table1D
     int begin = 1;
     int end = 0;
 
-    AMREX_GPU_HOST_DEVICE
-    constexpr Table1D () noexcept {}
+    constexpr Table1D () noexcept = default;
 
     template <class U=T, typename std::enable_if<std::is_const<U>::value,int>::type = 0>
     AMREX_GPU_HOST_DEVICE
@@ -81,8 +80,7 @@ struct Table2D
     GpuArray<int,2> begin{{1,1}};
     GpuArray<int,2> end{{0,0}};
 
-    AMREX_GPU_HOST_DEVICE
-    constexpr Table2D () noexcept {}
+    constexpr Table2D () noexcept = default;
 
     template <class U=T, typename std::enable_if<std::is_const<U>::value,int>::type = 0>
     AMREX_GPU_HOST_DEVICE
@@ -147,8 +145,7 @@ struct Table3D
     GpuArray<int,3> begin{{1,1,1}};
     GpuArray<int,3> end{{0,0,0}};
 
-    AMREX_GPU_HOST_DEVICE
-    constexpr Table3D () noexcept {}
+    constexpr Table3D () noexcept = default;
 
     template <class U=T, typename std::enable_if<std::is_const<U>::value,int>::type = 0>
     AMREX_GPU_HOST_DEVICE
@@ -219,8 +216,7 @@ struct Table4D
     GpuArray<int,4> begin{{1,1,1,1}};
     GpuArray<int,4> end{{0,0,0,0}};
 
-    AMREX_GPU_HOST_DEVICE
-    constexpr Table4D () noexcept {}
+    constexpr Table4D () noexcept = default;
 
     template <class U=T, typename std::enable_if<std::is_const<U>::value,int>::type = 0>
     AMREX_GPU_HOST_DEVICE
@@ -337,7 +333,7 @@ public:
                              std::conditional_t<N==3, Table3D<T const>,
                                                       Table4D<T const> > > >;
 
-    TableData () noexcept;
+    TableData () noexcept = default;
 
     explicit TableData (Arena* ar) noexcept;
 
@@ -351,11 +347,11 @@ public:
 
     ~TableData () noexcept;
 
-    constexpr int dim () const noexcept { return N; }
+    [[nodiscard]] constexpr int dim () const noexcept { return N; }
 
     void resize (Array<int,N> const& lo, Array<int,N> const& hi, Arena* ar = nullptr);
 
-    Long size () const noexcept;
+    [[nodiscard]] Long size () const noexcept;
 
     Array<int,N> const& lo () const noexcept { return m_lo; }
 
@@ -380,9 +376,6 @@ private:
     bool m_ptr_owner = false;
 };
 
-template <typename T, int N>
-TableData<T,N>::TableData () noexcept {}
-
 template <typename T, int N>
 TableData<T,N>::TableData (Array<int,N> const& lo, Array<int,N> const& hi, Arena* ar)
     : DataAllocator{ar}, m_lo(lo), m_hi(hi)
diff --git a/Src/LinearSolvers/AMReX_GMRES.H b/Src/LinearSolvers/AMReX_GMRES.H
new file mode 100644
index 00000000000..fd729c9e45d
--- /dev/null
+++ b/Src/LinearSolvers/AMReX_GMRES.H
@@ -0,0 +1,406 @@
+#ifndef AMREX_GMRES_H_
+#define AMREX_GMRES_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_BLProfiler.H>
+#include <AMReX_Print.H>
+#include <AMReX_TableData.H>
+#include <AMReX_Vector.H>
+#include <cmath>
+#include <limits>
+#include <memory>
+
+namespace amrex {
+
+/**
+ * \brief GMRES
+ *
+ * This class implements the GMRES algorithm. The template parameter V is
+ * for a linear algebra vector class. For example, it could be
+ * amrex::MultiFab. The other template parameter M is for a linear operator
+ * class with a number of required member functions. Note that conceptually
+ * M contains a matrix. However, it does not mean it needs to have a member
+ * variable storing the matrix, because GMRES only needs the matrix vector
+ * product, not the matrix itself.
+ *
+ * \tparam V linear algebra vector. It must be default constructible, move
+ *           constructible, and move assignable.
+ * \tparam M linear operator. A list of required member functions for M is
+ *           shown below. Here RT (typename V::value_type) is either double
+ *           or float.
+ *             - void apply(V& lhs, V const& rhs)\n
+ *               lhs = L(rhs), where L is the linear operator performing matrix
+ *               vector product.
+ *             - void assign(V& lhs, V const& rhs)\n
+ *               lhs = rhs.
+ *             - RT dotProduct(V const& v1, V const& v2)\n
+ *               returns v1 * v2.
+ *             - void increment(V& lhs, V const& rhs, RT a)\n
+ *               lhs += a * rhs.
+ *             - void linComb(V& lhs, RT a, V const& rhs_a, RT b, V const& rhs_b)\n
+ *               lhs = a * rhs_a + b * rhs_b.
+ *             - V makeVecRHS()\n
+ *               returns a V object that is suitable as RHS in M x = b. The reason
+ *               we distinguish between LHS and RHS is M might need the distinction
+ *               for efficiency. For example, if V is MultiFab, we might need the x
+ *               in the LHS of M x = b to have ghost cells for efficiency, whereas
+ *               no ghost cells are needed for the RHS (i.e., b).
+ *             - V makeVecLHS()\n
+ *               returns a V object that is suitable as LHS in M x = b. See the
+ *               description for makeVecRHS for more details.
+ *             - RT norm2(V const& v)\n
+ *               returns the 2-norm of v.
+ *             - void precond(V& lhs, V const& rhs)\n
+ *               applies preconditioner to rhs. If there is no preconditioner,
+ *               this function should do lhs = rhs.
+ *             - void setVal(V& v, RT value)\n
+ *               v = value.
+ */
+template <typename V, typename M>
+class GMRES
+{
+public:
+
+    using RT = typename V::value_type; // double or float
+
+    GMRES ();
+
+    //! Defines with a reference to M. It's the user's responsibility to
+    //! keep the M object alive for GMRES to be functional. This function
+    //! must be called before solve() can be called.
+    void define (M& linop);
+
+    /**
+     * \brief Solve the linear system
+     *
+     * \param a_sol     unknowns, i.e., x in A x = b.
+     * \param a_rhs     RHS, i.e., b in A x = b.
+     * \param a_tol_rel relative tolerance.
+     * \param a_tol_abs absolute tolerance.
+     * \patam a_its     optional argument specifying the maximum number of iterations.
+     */
+    void solve (V& a_sol, V const& a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its=-1);
+
+    //! Sets verbosity.
+    void setVerbose (int v) { m_verbose = v; }
+
+    //! Sets restart length. The default is 30.
+    void setRestartLength (int rl);
+
+    //! Gets the number of iterations.
+    [[nodiscard]] int getNumIters () const { return m_its; }
+
+    //! Gets the solver status.
+    [[nodiscard]] int getStatus () const { return m_status; }
+
+    //! Gets the 2-norm of the residual.
+    [[nodiscard]] RT getResidualNorm () const { return m_res; }
+
+private:
+    void clear ();
+    void allocate_scratch ();
+    void cycle (V& a_xx, int& a_status, int& a_itcount, RT& a_rnorm0);
+    void build_solution (V& a_xx, int it);
+    void compute_residual (V& a_rr, V const& a_xx, V const& a_bb);
+
+    bool converged (RT r0, RT r) const;
+
+    void gram_schmidt_orthogonalization (int it);
+    void update_hessenberg (int it, bool happyend, RT& res);
+
+    int m_verbose = 0;
+    int m_maxiter = 2000;
+    int m_its = 0;
+    int m_status = -1;
+    int m_restrtlen = 30;
+    RT m_res = std::numeric_limits<RT>::max();
+    RT m_rtol = RT(0);
+    RT m_atol = RT(0);
+    Vector<RT> m_hh_1d;
+    Vector<RT> m_hes_1d;
+    Table2D<RT> m_hh;
+    Table2D<RT> m_hes;
+    Vector<RT> m_grs;
+    Vector<RT> m_cc;
+    Vector<RT> m_ss;
+    std::unique_ptr<V> m_v_tmp_rhs;
+    std::unique_ptr<V> m_v_tmp_lhs;
+    Vector<V> m_vv;
+    M* m_linop = nullptr;
+};
+
+template <typename V, typename M>
+GMRES<V,M>::GMRES ()
+{
+    allocate_scratch();
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::allocate_scratch ()
+{
+    int rs = m_restrtlen;
+
+    m_hh_1d.resize(std::size_t(rs + 2) * (rs + 1));
+    m_hh = Table2D<RT>(m_hh_1d.data(), {0,0}, {rs+1,rs}); // (0:rs+1,0:rs)
+
+    m_hes_1d.resize(std::size_t(rs + 2) * (rs + 1));
+    m_hes = Table2D<RT>(m_hes_1d.data(), {0,0}, {rs+1,rs}); // (0:rs+1,0:rs)
+
+    m_grs.resize(rs + 2);
+    m_cc.resize(rs + 1);
+    m_ss.resize(rs + 1);
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::setRestartLength (int rl)
+{
+    if (m_restrtlen != rl) {
+        m_restrtlen = rl;
+        allocate_scratch();
+        m_vv.clear();
+    }
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::define (M& linop)
+{
+    clear();
+    m_linop = &linop;
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::clear ()
+{
+    m_its = 0;
+    m_status = -1;
+    m_res = std::numeric_limits<RT>::max();
+    m_v_tmp_rhs.reset();
+    m_v_tmp_lhs.reset();
+    m_vv.clear();
+    m_linop = nullptr;
+}
+
+template <typename V, typename M>
+bool GMRES<V,M>::converged (RT r0, RT r) const
+{
+    return (r < r0*m_rtol) || (r < m_atol);
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::solve (V& a_sol, V const& a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its)
+{
+    BL_PROFILE("GMRES::solve()");
+
+    AMREX_ALWAYS_ASSERT(m_linop != nullptr);
+
+    auto t0 = amrex::second();
+
+    if (m_v_tmp_rhs == nullptr) {
+        m_v_tmp_rhs = std::make_unique<V>(m_linop->makeVecRHS());
+    }
+    if (m_v_tmp_lhs == nullptr) {
+        m_v_tmp_lhs = std::make_unique<V>(m_linop->makeVecLHS());
+    }
+    if (m_vv.empty()) {
+        m_vv.resize(m_restrtlen+1);
+        for (auto& v : m_vv) {
+            v = m_linop->makeVecRHS();
+        }
+    }
+
+    m_rtol = a_tol_rel;
+    m_atol = a_tol_abs;
+
+    if (a_its < 0) { a_its = m_maxiter; }
+
+    auto rnorm0 = RT(0);
+
+    m_linop->assign(m_vv[0], a_rhs);
+    m_linop->setVal(a_sol, RT(0.0));
+
+    m_its = 0;
+    m_status = -1;
+    cycle(a_sol, m_status, m_its, rnorm0);
+
+    while (m_status == -1 && m_its < a_its) {
+        compute_residual(m_vv[0], a_sol, a_rhs);
+        cycle(a_sol, m_status, m_its, rnorm0);
+    }
+
+    if (m_status == -1 && m_its >= a_its) { m_status = 1; }
+
+    m_v_tmp_rhs.reset();
+    m_v_tmp_lhs.reset();
+    m_vv.clear();
+
+    auto t1 = amrex::second();
+    if (m_verbose > 0) {
+        amrex::Print() << "GMRES: Solve Time = " << t1-t0 << std::endl;
+    }
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::cycle (V& a_xx, int& a_status, int& a_itcount, RT& a_rnorm0)
+{
+    BL_PROFILE("GMREA::cycle()");
+
+    m_res = m_linop->norm2(m_vv[0]);
+    m_grs[0] = m_res;
+
+    if (m_res == RT(0.0)) {
+        a_status = 0;
+        return;
+    }
+
+    m_linop->scale(m_vv[0], RT(1.0)/m_res);
+
+    if (a_itcount == 0) { a_rnorm0 = m_res; }
+
+    a_status = converged(a_rnorm0,m_res) ? 0 : -1;
+
+    int it = 0;
+    while (it < m_restrtlen && a_itcount < m_maxiter)
+    {
+        if (m_verbose > 1) {
+            amrex::Print() << "GMRES: iter = " << a_itcount
+                           << ", residual = " << m_res << ", " << m_res/a_rnorm0
+                           << " (rel.)\n";
+        }
+
+        if (a_status == 0) { break; }
+
+        auto const& vv_it  = m_vv[it  ];
+        auto      & vv_it1 = m_vv[it+1];
+
+        m_linop->precond(*m_v_tmp_lhs, vv_it);
+        m_linop->apply(vv_it1, *m_v_tmp_lhs);
+
+        gram_schmidt_orthogonalization(it);
+
+        auto tt = m_linop->norm2(vv_it1);
+
+        auto const small = RT((sizeof(RT) == 8) ? 1.e-99 : 1.e-30);
+        bool happyend = (tt < small);
+        if (!happyend) {
+            m_linop->scale(vv_it1, RT(1.0)/tt);
+        }
+
+        m_hh (it+1,it) = tt;
+        m_hes(it+1,it) = tt;
+
+        update_hessenberg(it, happyend, m_res);
+
+        ++it;
+        ++a_itcount;
+        a_status = converged(a_rnorm0, m_res) ? 0 : -1;
+        if (happyend) { break; }
+    }
+
+    if ((m_verbose > 1) && (a_status != 0 || a_itcount >= m_maxiter)) {
+        amrex::Print() << "GMRES: iter = " << a_itcount
+                       << ", residual = " << m_res << ", " << m_res/a_rnorm0
+                       << " (rel.)\n";
+    }
+
+    build_solution(a_xx, it-1);
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::gram_schmidt_orthogonalization (int const it)
+{
+    // Two unmodified Gram-Schmidt Orthogonalization
+
+    BL_PROFILE("GMRES::GramSchmidt");
+
+    auto& vv_1 = m_vv[it+1];
+
+    Vector<RT> lhh(it+1);
+
+    for (int j = 0; j <= it; ++j) {
+        m_hh (j,it) = RT(0.0);
+        m_hes(j,it) = RT(0.0);
+    }
+
+    for (int ncnt = 0; ncnt < 2 ; ++ncnt)
+    {
+        for (int j = 0; j <= it; ++j) {
+            lhh[j] = m_linop->dotProduct(vv_1, m_vv[j]);
+        }
+
+        for (int j = 0; j <= it; ++j) {
+            m_linop->increment(vv_1, m_vv[j], -lhh[j]);
+            m_hh (j,it) += lhh[j];
+            m_hes(j,it) -= lhh[j];
+        }
+    }
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::update_hessenberg (int it, bool happyend, RT& res)
+{
+    BL_PROFILE("GMRES::update_hessenberg()");
+
+    for (int j = 1; j <= it; ++j) {
+        auto tt = m_hh(j-1,it);
+        m_hh(j-1,it) = m_cc[j-1] * tt         + m_ss[j-1] * m_hh(j,it);
+        m_hh(j  ,it) = m_cc[j-1] * m_hh(j,it) - m_ss[j-1] * tt;
+    }
+
+    if (!happyend)
+    {
+        auto tt = std::sqrt(m_hh(it,it)*m_hh(it,it) + m_hh(it+1,it)*m_hh(it+1,it));
+        m_cc[it] = m_hh(it  ,it) / tt;
+        m_ss[it] = m_hh(it+1,it) / tt;
+        m_grs[it+1] = - (m_ss[it] * m_grs[it]);
+        m_grs[it  ] =    m_cc[it] * m_grs[it];
+        m_hh(it,it) = m_cc[it] * m_hh(it,it) + m_ss[it] * m_hh(it+1,it);
+        res = std::abs(m_grs[it+1]);
+    }
+    else
+    {
+        res = RT(0.0);
+    }
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::build_solution (V& a_xx, int const it)
+{
+    BL_PROFILE("GMRES:build_solution()");
+
+    if (it < 0) { return; }
+
+    if (m_hh(it,it) != RT(0.0)) {
+        m_grs[it] /= m_hh(it,it);
+    } else {
+        m_grs[it] = RT(0.0);
+    }
+
+    for (int ii = 1; ii <= it; ++ii) {
+        int k = it - ii;
+        auto tt = m_grs[k];
+        for (int j = k+1; j <= it; ++j) {
+            tt -= m_hh(k,j) * m_grs[j];
+        }
+        m_grs[k] = tt / m_hh(k,k);
+    }
+
+    m_linop->setVal(*m_v_tmp_rhs, RT(0.0));
+    for (int ii = 0; ii < it+1; ++ii) {
+        m_linop->increment(*m_v_tmp_rhs, m_vv[ii], m_grs[ii]);
+    }
+
+    m_linop->precond(*m_v_tmp_lhs, *m_v_tmp_rhs);
+    m_linop->increment(a_xx, *m_v_tmp_lhs, RT(1.0));
+}
+
+template <typename V, typename M>
+void GMRES<V,M>::compute_residual (V& a_rr, V const& a_xx, V const& a_bb)
+{
+    BL_PROFILE("GMRES::compute_residual()");
+    m_linop->assign(*m_v_tmp_lhs, a_xx);
+    m_linop->apply(*m_v_tmp_rhs, *m_v_tmp_lhs);
+    m_linop->linComb(a_rr, RT(1.0), a_bb, RT(-1.0), *m_v_tmp_rhs);
+}
+
+}
+#endif
diff --git a/Src/LinearSolvers/AMReX_GMRES_MLMG.H b/Src/LinearSolvers/AMReX_GMRES_MLMG.H
new file mode 100644
index 00000000000..5106afde378
--- /dev/null
+++ b/Src/LinearSolvers/AMReX_GMRES_MLMG.H
@@ -0,0 +1,147 @@
+#ifndef AMREX_GMRES_MLMG_H_
+#define AMREX_GMRES_MLMG_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_MLMG.H>
+#include <utility>
+
+namespace amrex {
+
+//! Wrapping MLMG as a matrix operator for GMRES
+template <typename M>
+class GMRESMLMGT
+{
+public:
+    using MF = typename M::MFType; // typically MultiFab
+    using RT = typename MF::value_type; // double or float
+
+    explicit GMRESMLMGT (M& mlmg);
+
+    //! Make MultiFab without ghost cells
+    MF makeVecRHS () const;
+
+    //! Make MultiFab with ghost cells and set ghost cells to zero
+    MF makeVecLHS () const;
+
+    RT norm2 (MF const& mf) const;
+
+    static void scale (MF& mf, RT scale_factor);
+
+    RT dotProduct (MF const& mf1, MF const& mf2) const;
+
+    //! lhs = value
+    static void setVal (MF& lhs, RT value);
+
+    //! lhs = rhs
+    static void assign (MF& lhs, MF const& rhs);
+
+    //! lhs += a*rhs
+    static void increment (MF& lhs, MF const& rhs, RT a);
+
+    //! lhs = a*rhs_a + b*rhs_b
+    static void linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b);
+
+    //! lhs = L(rhs)
+    void apply (MF& lhs, MF const& rhs) const;
+
+    void precond (MF& lhs, MF const& rhs) const;
+
+    bool usePrecond (bool new_flag) { return std::exchange(m_use_precond, new_flag); }
+
+private:
+    M& m_mlmg;
+    MLLinOpT<MF>& m_linop;
+    bool m_use_precond = false;
+};
+
+template <typename M>
+GMRESMLMGT<M>::GMRESMLMGT (M& mlmg)
+    : m_mlmg(mlmg), m_linop(mlmg.getLinOp())
+{
+    m_mlmg.prepareLinOp();
+}
+
+template <typename M>
+auto GMRESMLMGT<M>::makeVecRHS () const -> MF
+{
+    return m_linop.make(0, 0, IntVect(0));
+}
+
+template <typename M>
+auto GMRESMLMGT<M>::makeVecLHS () const -> MF
+{
+    auto mf = m_linop.make(0, 0, IntVect(1));
+    mf.setBndry(0);
+    return mf;
+}
+
+template <typename M>
+auto GMRESMLMGT<M>::norm2 (MF const& mf) const -> RT
+{
+    auto r = m_linop.xdoty(0, 0, mf, mf, false);
+    return std::sqrt(r);
+}
+
+template <typename M>
+void GMRESMLMGT<M>::scale (MF& mf, RT scale_factor)
+{
+    mf.mult(scale_factor, 0, mf.nComp());
+}
+
+template <typename M>
+auto GMRESMLMGT<M>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
+{
+    return m_linop.xdoty(0, 0, mf1, mf2, false);
+}
+
+template <typename M>
+void GMRESMLMGT<M>::setVal (MF& lhs, RT value)
+{
+    lhs.setVal(value);
+}
+
+template <typename M>
+void GMRESMLMGT<M>::assign (MF& lhs, MF const& rhs)
+{
+    MF::Copy(lhs, rhs, 0, 0, lhs.nComp(), IntVect(0));
+}
+
+template <typename M>
+void GMRESMLMGT<M>::increment (MF& lhs, MF const& rhs, RT a)
+{
+    MF::Saxpy(lhs, a, rhs, 0, 0, lhs.nComp(), IntVect(0));
+}
+
+template <typename M>
+void GMRESMLMGT<M>::linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const& rhs_b)
+{
+    MF::LinComb(lhs, a, rhs_a, 0, b, rhs_b, 0, 0, lhs.nComp(), IntVect(0));
+}
+
+template <typename M>
+void GMRESMLMGT<M>::apply (MF& lhs, MF const& rhs) const
+{
+    m_linop.apply(0, 0, lhs, const_cast<MF&>(rhs),
+                  MLLinOpT<MF>::BCMode::Homogeneous,
+                  MLLinOpT<MF>::StateMode::Correction);
+}
+
+template <typename M>
+void GMRESMLMGT<M>::precond (MF& lhs, MF const& rhs) const
+{
+    if (m_use_precond) {
+        // for now, let's just do some smoothing
+        lhs.setVal(RT(0.0));
+        for (int m = 0; m < 4; ++m) {
+            m_linop.smooth(0, 0, lhs, rhs, (m==0) ? true : false);
+        }
+    } else {
+        amrex::Copy(lhs, rhs, 0, 0, lhs.nComp(), IntVect(0));
+    }
+}
+
+using GMRESMLMG = GMRESMLMGT<MLMG>;
+
+}
+
+#endif
diff --git a/Src/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt
index c2851d49959..76f75e06123 100644
--- a/Src/LinearSolvers/CMakeLists.txt
+++ b/Src/LinearSolvers/CMakeLists.txt
@@ -2,6 +2,7 @@ foreach(D IN LISTS AMReX_SPACEDIM)
     #
     # Sources in subdirectory MLMG
     #
+    target_include_directories(amrex_${D}d PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>)
     target_include_directories(amrex_${D}d PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/MLMG>)
 
     target_sources(amrex_${D}d
@@ -53,6 +54,8 @@ foreach(D IN LISTS AMReX_SPACEDIM)
        MLMG/AMReX_MLNodeABecLaplacian.cpp
        MLMG/AMReX_MLNodeABecLap_K.H
        MLMG/AMReX_MLNodeABecLap_${D}D_K.H
+       AMReX_GMRES.H
+       AMReX_GMRES_MLMG.H
        )
 
     if (D EQUAL 3)
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H
index f0dca07f3ab..5fc8de10022 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H
@@ -98,6 +98,7 @@ template <typename T> class MLMGT;
 template <typename T> class MLCGSolverT;
 template <typename T> class MLPoissonT;
 template <typename T> class MLABecLaplacianT;
+template <typename T> class GMRESMLMGT;
 
 template <typename MF>
 class MLLinOpT
@@ -108,7 +109,9 @@ public:
     template <typename T> friend class MLCGSolverT;
     template <typename T> friend class MLPoissonT;
     template <typename T> friend class MLABecLaplacianT;
+    template <typename T> friend class GMRESMLMGT;
 
+    using MFType = MF;
     using FAB = typename LinOpData<MF>::fab_type;
     using RT  = typename LinOpData<MF>::value_type;
 
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H
index 9bfc2f0007d..46069b7a26b 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H
@@ -21,6 +21,7 @@ public:
 
     template <typename T> friend class MLCGSolverT;
 
+    using MFType = MF;
     using FAB = typename MLLinOpT<MF>::FAB;
     using RT  = typename MLLinOpT<MF>::RT;
 
@@ -167,6 +168,8 @@ public:
 
     void prepareForNSolve ();
 
+    void prepareLinOp ();
+
     void oneIter (int iter);
 
     void miniCycle (int amrlev);
@@ -217,6 +220,8 @@ public:
     [[nodiscard]] int getNumIters () const noexcept { return m_iter_fine_resnorm0.size(); }
     [[nodiscard]] Vector<int> const& getNumCGIters () const noexcept { return m_niters_cg; }
 
+    MLLinOpT<MF>& getLinOp () { return linop; }
+
 private:
 
     bool throw_exception = false;
@@ -794,12 +799,7 @@ MLMGT<MF>::compResidual (const Vector<MF*>& a_res, const Vector<MF*>& a_sol,
         }
     }
 
-    if (!linop_prepared) {
-        linop.prepareForSolve();
-        linop_prepared = true;
-    } else if (linop.needsUpdate()) {
-        linop.update();
-    }
+    prepareLinOp();
 
     const auto& amrrr = linop.AMRRefRatio();
 
@@ -880,12 +880,7 @@ MLMGT<MF>::apply (const Vector<MF*>& out, const Vector<MF*>& a_in)
         setVal(rh[alev], RT(0.0));
     }
 
-    if (!linop_prepared) {
-        linop.prepareForSolve();
-        linop_prepared = true;
-    } else if (linop.needsUpdate()) {
-        linop.update();
-    }
+    prepareLinOp();
 
     for (int alev = 0; alev < namrlevs; ++alev) {
         linop.applyInhomogNeumannTerm(alev, rh[alev]);
@@ -1106,6 +1101,18 @@ MLMGT<MF>::prepareForSolve (Vector<AMF*> const& a_sol, Vector<AMF const*> const&
     }
 }
 
+template <typename MF>
+void
+MLMGT<MF>::prepareLinOp ()
+{
+    if (!linop_prepared) {
+        linop.prepareForSolve();
+        linop_prepared = true;
+    } else if (linop.needsUpdate()) {
+        linop.update();
+    }
+}
+
 template <typename MF>
 void
 MLMGT<MF>::prepareForNSolve ()
diff --git a/Src/LinearSolvers/MLMG/Make.package b/Src/LinearSolvers/MLMG/Make.package
index d66d64ec0eb..769ccad4d26 100644
--- a/Src/LinearSolvers/MLMG/Make.package
+++ b/Src/LinearSolvers/MLMG/Make.package
@@ -1,3 +1,5 @@
+ifndef AMREX_MLMG_MAKE
+       AMREX_MLMG_MAKE := 1
 
 CEXE_sources += AMReX_MLMG.cpp
 
@@ -87,3 +89,5 @@ endif
 
 VPATH_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/MLMG
 INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/MLMG
+
+endif
diff --git a/Src/LinearSolvers/Make.package b/Src/LinearSolvers/Make.package
new file mode 100644
index 00000000000..c818714cbfd
--- /dev/null
+++ b/Src/LinearSolvers/Make.package
@@ -0,0 +1,9 @@
+CEXE_headers += AMReX_GMRES.H AMReX_GMRES_MLMG.H
+
+VPATH_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers
+INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers
+
+include $(AMREX_HOME)/Src/LinearSolvers/MLMG/Make.package
+ifeq ($(DIM),3)
+  include $(AMREX_HOME)/Src/LinearSolvers/OpenBC/Make.package
+endif
diff --git a/Src/LinearSolvers/OpenBC/Make.package b/Src/LinearSolvers/OpenBC/Make.package
index 5fc39f69371..b2f83ed61e8 100644
--- a/Src/LinearSolvers/OpenBC/Make.package
+++ b/Src/LinearSolvers/OpenBC/Make.package
@@ -1,6 +1,10 @@
+ifndef AMREX_OPENBC_MAKE
+       AMREX_OPENBC_MAKE := 1
 
 CEXE_headers += AMReX_OpenBC.H AMReX_OpenBC_K.H
 CEXE_sources += AMReX_OpenBC.cpp
 
 VPATH_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/OpenBC
 INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/OpenBC
+
+endif
diff --git a/Tests/LinearSolvers/ABecLaplacian_C/GNUmakefile b/Tests/LinearSolvers/ABecLaplacian_C/GNUmakefile
index 329fb9afcf2..e26200bd94f 100644
--- a/Tests/LinearSolvers/ABecLaplacian_C/GNUmakefile
+++ b/Tests/LinearSolvers/ABecLaplacian_C/GNUmakefile
@@ -16,7 +16,7 @@ include $(AMREX_HOME)/Tools/GNUMake/Make.defs
 
 include ./Make.package
 
-Pdirs 	:= Base Boundary LinearSolvers/MLMG
+Pdirs 	:= Base Boundary LinearSolvers
 
 Ppack	+= $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package)
 
diff --git a/Tests/LinearSolvers/ABecLaplacian_C/MyTest.H b/Tests/LinearSolvers/ABecLaplacian_C/MyTest.H
index c1ed7ba4c3d..97c928a88da 100644
--- a/Tests/LinearSolvers/ABecLaplacian_C/MyTest.H
+++ b/Tests/LinearSolvers/ABecLaplacian_C/MyTest.H
@@ -30,6 +30,7 @@ private:
     void solveABecLaplacian ();
     void solveABecLaplacianInhomNeumann ();
     void solveNodeABecLaplacian ();
+    void solveABecLaplacianGMRES ();
 
     int max_level = 1;
     int ref_ratio = 2;
@@ -56,6 +57,9 @@ private:
     bool use_hypre = false;
     bool use_petsc = false;
 
+    // GMRES
+    bool use_gmres = false;
+
 #ifdef AMREX_USE_HYPRE
     int hypre_interface_i = 1;  // 1. structed, 2. semi-structed, 3. ij
     amrex::Hypre::Interface hypre_interface = amrex::Hypre::Interface::structed;
diff --git a/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp b/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp
index 9900ce43d7f..84fe99cb9cb 100644
--- a/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp
+++ b/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp
@@ -1,5 +1,7 @@
 #include "MyTest.H"
 
+#include <AMReX_GMRES.H>
+#include <AMReX_GMRES_MLMG.H>
 #include <AMReX_MLNodeABecLaplacian.H>
 #include <AMReX_MLABecLaplacian.H>
 #include <AMReX_MLPoisson.H>
@@ -20,7 +22,11 @@ MyTest::solve ()
     if (prob_type == 1) {
         solvePoisson();
     } else if (prob_type == 2) {
-        solveABecLaplacian();
+        if (use_gmres) {
+            solveABecLaplacianGMRES();
+        } else {
+            solveABecLaplacian();
+        }
     } else if (prob_type == 3) {
         solveABecLaplacianInhomNeumann();
     } else if (prob_type == 4) {
@@ -460,6 +466,86 @@ MyTest::solveNodeABecLaplacian ()
     }
 }
 
+void
+MyTest::solveABecLaplacianGMRES ()
+{
+    LPInfo info;
+    info.setMaxCoarseningLevel(0);
+
+    const auto tol_rel = Real(1.e-6);
+    const auto tol_abs = Real(0.0);
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(composite_solve == false,
+       "solveABecLaplacianGMRES does not support composite solve");
+
+    const auto nlevels = static_cast<int>(geom.size());
+
+    for (int ilev = 0; ilev < nlevels; ++ilev)
+    {
+        MLABecLaplacian mlabec({geom[ilev]}, {grids[ilev]}, {dmap[ilev]}, info);
+
+        mlabec.setMaxOrder(linop_maxorder);
+
+        // This is a 3d problem with homogeneous Neumann BC
+        mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann,
+                                         LinOpBCType::Neumann,
+                                         LinOpBCType::Neumann)},
+            {AMREX_D_DECL(LinOpBCType::Neumann,
+                          LinOpBCType::Neumann,
+                          LinOpBCType::Neumann)});
+
+        if (ilev > 0) {
+            mlabec.setCoarseFineBC(&solution[ilev-1], ref_ratio);
+        }
+
+        // for problem with pure homogeneous Neumann BC, we could pass a nullptr
+        mlabec.setLevelBC(0, nullptr);
+
+        mlabec.setScalars(ascalar, bscalar);
+
+        mlabec.setACoeffs(0, acoef[ilev]);
+
+        Array<MultiFab,AMREX_SPACEDIM> face_bcoef;
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
+        {
+            const BoxArray& ba = amrex::convert(bcoef[ilev].boxArray(),
+                                                IntVect::TheDimensionVector(idim));
+            face_bcoef[idim].define(ba, bcoef[ilev].DistributionMap(), 1, 0);
+        }
+        amrex::average_cellcenter_to_face(GetArrOfPtrs(face_bcoef),
+                                          bcoef[ilev], geom[ilev]);
+        mlabec.setBCoeffs(0, amrex::GetArrOfConstPtrs(face_bcoef));
+
+        MultiFab res(rhs[ilev].boxArray(), rhs[ilev].DistributionMap(), 1, 0);
+
+        MLMG mlmg(mlabec);
+        mlmg.setVerbose(verbose);
+        mlmg.apply({&res}, {&solution[ilev]}); // res = L(sol)
+
+        MultiFab::Subtract(res, rhs[ilev], 0, 0, 1, 0); // now res = L(sol) - rhs
+
+        MultiFab cor(rhs[ilev].boxArray(), rhs[ilev].DistributionMap(), 1, 0);
+
+        using M = GMRESMLMG;
+        M mat(mlmg);
+        mat.usePrecond(true);
+
+        GMRES<MultiFab,M> gmres;
+        gmres.setVerbose(verbose);
+        gmres.define(mat);
+        gmres.solve(cor, res, tol_rel, tol_abs); // solve L(cor) = res
+
+        MultiFab::Subtract(solution[ilev], cor, 0, 0, 1, 0);
+
+        mlmg.apply({&res}, {&solution[ilev]}); // res = L(sol)
+        MultiFab::Subtract(res, rhs[ilev], 0, 0, 1, 0); // now res = L(sol) - rhs
+        if (verbose) {
+            amrex::Print() << "Final residual = " << res.norminf(0)
+                           << " " << res.norm1(0) << " " << res.norm2(0) << std::endl;
+        }
+    }
+}
+
 void
 MyTest::readParameters ()
 {
@@ -484,6 +570,9 @@ MyTest::readParameters ()
     pp.query("max_coarsening_level", max_coarsening_level);
     pp.query("max_semicoarsening_level", max_semicoarsening_level);
 
+    pp.query("use_gmres", use_gmres);
+    AMREX_ALWAYS_ASSERT(use_gmres == false || prob_type == 2);
+
 #ifdef AMREX_USE_HYPRE
     pp.query("use_hypre", use_hypre);
     pp.query("hypre_interface", hypre_interface_i);
diff --git a/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp b/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp
index 4473f978a85..6dbb2c55b59 100644
--- a/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp
+++ b/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp
@@ -77,7 +77,8 @@ MyTest::writePlotfile () const
             Real dvol = AMREX_D_TERM(dx[0],*dx[1],*dx[2]);
             amrex::Print() << "Level " << ilev
                            << " max-norm error: " << plotmf[ilev].norminf(3)
-                           << " 1-norm error: " << plotmf[ilev].norm1(3)*dvol << std::endl;
+                           << " 1-norm error: " << plotmf[ilev].norm1(3)*dvol
+                           << " 2-norm error: " << plotmf[ilev].norm2(3)*dvol << std::endl;
         }
 
         WriteMultiLevelPlotfile("plot", nlevels, amrex::GetVecOfConstPtrs(plotmf),
diff --git a/Tests/LinearSolvers/ABecLaplacian_C/inputs.gmres b/Tests/LinearSolvers/ABecLaplacian_C/inputs.gmres
new file mode 100644
index 00000000000..17e7767dffa
--- /dev/null
+++ b/Tests/LinearSolvers/ABecLaplacian_C/inputs.gmres
@@ -0,0 +1,17 @@
+
+max_level = 0
+ref_ratio = 2
+n_cell = 128
+max_grid_size = 64
+
+composite_solve = 0   # composite solve or level by level?
+
+prob_type = 2
+
+use_gmres = 1
+
+verbose = 2
+
+amrex.fpe_trap_invalid = 1
+amrex.fpe_trap_zero = 1
+amrex.fpe_trap_overflow = 1