diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
index 974ad14d..bd72c073 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/AMReX_GMRES_Poisson.H
@@ -21,7 +21,7 @@ public:
     using RT = amrex::Real; // typename MF::RT; // double or float
     using GM = GMRES<MF,GMRESPOISSONT<MF>>;
 
-    explicit GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm);
+    explicit GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom);
 
     /**
      * \brief Solve the linear system
@@ -36,30 +36,9 @@ public:
     //! Sets verbosity.
     void setVerbose (int v) { m_gmres.setVerbose(v); }
 
-    //! Sets the max number of iterations
-    void setMaxIters (int niters) { m_gmres.setMaxIters(niters); }
-
-    //! Gets the number of iterations.
-    [[nodiscard]] int getNumIters () const { return m_gmres.getNumIters(); }
-
-    //! Gets the 2-norm of the residual.
-    [[nodiscard]] RT getResidualNorm () const { return m_gmres.getResidualNorm(); }
-
     //! Get the GMRES object.
     GM& getGMRES () { return m_gmres; }
 
-    /**
-     * \brief Set MLMG's multiplicative property of zero
-     *
-     * This should NOT be called unless MLMG has the multiplicative property
-     * of zero. MLMG is not a matrix, and usually does not have the
-     * properties of a matrix such as the multiplicative property of zero
-     * (i.e., M*0=0) because how domain boundary conditions are
-     * handled. However, if MLMG has the property of zero, calling this
-     * function with true can have a small performance benefit.
-     */
-    void setPropertyOfZero (bool b) { m_prop_zero = b; }
-
     //! Make MultiFab without ghost cells
     MF makeVecRHS () const;
 
@@ -85,28 +64,24 @@ public:
     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 apply (MF& lhs, MF& rhs) const;
 
     void precond (MF& lhs, MF const& rhs) const;
 
     //! Control whether or not to use MLMG as preconditioner.
     bool usePrecond (bool new_flag) { return std::exchange(m_use_precond, new_flag); }
 
-    //! Set the number of MLMG preconditioner iterations per GMRES iteration.
-    void setPrecondNumIters (int precond_niters) { m_precond_niters = precond_niters; }
-
 private:
     GM m_gmres;
     BoxArray m_ba;
     DistributionMapping m_dm;
-    bool m_use_precond = true;
-    bool m_prop_zero = false;
-    int m_precond_niters = 1;
+    Geometry m_geom;
+    bool m_use_precond = false;
 };
 
 template <typename MF>
-GMRESPOISSONT<MF>::GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm)
-    : m_ba(ba), m_dm(dm)
+GMRESPOISSONT<MF>::GMRESPOISSONT (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom)
+    : m_ba(ba), m_dm(dm), m_geom(geom)
 {
     m_gmres.define(*this);
 }
@@ -128,8 +103,8 @@ auto GMRESPOISSONT<MF>::makeVecLHS () const -> MF
 template <typename MF>
 auto GMRESPOISSONT<MF>::norm2 (MF const& mf) const -> RT
 {
-//    auto r = m_linop->xdoty(0, 0, mf, mf, false);
-//    return std::sqrt(r);
+    auto r = MultiFab::Dot(mf,0,1,0);
+    return std::sqrt(r);
     return 0.;
 }
 
@@ -142,8 +117,7 @@ void GMRESPOISSONT<MF>::scale (MF& mf, RT scale_factor)
 template <typename MF>
 auto GMRESPOISSONT<MF>::dotProduct (MF const& mf1, MF const& mf2) const -> RT
 {
-//    return m_linop->xdoty(0, 0, mf1, mf2, false);
-    return 0.;
+    return MultiFab::Dot(mf1,0,mf2,0,1,0);
 }
 
 template <typename MF>
@@ -171,11 +145,33 @@ void GMRESPOISSONT<MF>::linComb (MF& lhs, RT a, MF const& rhs_a, RT b, MF const&
 }
 
 template <typename MF>
-void GMRESPOISSONT<MF>::apply (MF& lhs, MF const& rhs) const
+void GMRESPOISSONT<MF>::apply (MF& lhs, MF& rhs) const
 {
-//    m_linop->apply(0, 0, lhs, const_cast<MF&>(rhs),
-//                   MLLinOpT<MF>::BCMode::Homogeneous,
-//                   MLLinOpT<MF>::StateMode::Correction);
+    // apply matrix to rhs for output lhs
+    rhs.FillBoundary(m_geom.periodicity());
+    
+    const GpuArray<Real, AMREX_SPACEDIM> dx = m_geom.CellSizeArray();
+
+    for ( MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        
+        const Box& bx = mfi.tilebox();
+
+        const Array4<const Real> & rhs_p = rhs.array(mfi);
+        const Array4<      Real> & lhs_p = lhs.array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            lhs_p(i,j,k) = ( rhs_p(i+1,j,k) - 2.*rhs_p(i,j,k) + rhs_p(i-1,j,k) ) / (dx[0]*dx[0])
+                         + ( rhs_p(i,j+1,k) - 2.*rhs_p(i,j,k) + rhs_p(i,j-1,k) ) / (dx[1]*dx[1])
+#if (AMREX_SPACEDIM == 3)
+                         + ( rhs_p(i,j,k+1) - 2.*rhs_p(i,j,k) + rhs_p(i,j,k-1) ) / (dx[2]*dx[2])
+#endif
+                ;
+        });
+
+    }
+
+    
 }
 
 template <typename MF>
@@ -183,7 +179,6 @@ void GMRESPOISSONT<MF>::precond (MF& lhs, MF const& rhs) const
 {
     if (m_use_precond) {
 
-
     } else {
         LocalCopy(lhs, rhs, 0, 0, nComp(lhs), IntVect(0));
     }
@@ -201,6 +196,12 @@ void GMRESPOISSONT<MF>::solve (MF& a_sol, MF const& a_rhs, RT a_tol_rel, RT a_to
     m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
     increment(a_sol, cor, RT(-1)); // sol = sol - cor
     */
+    auto res = makeVecRHS();
+    apply(res,a_sol);
+    increment(res, a_rhs, RT(-1)); // res = L(sol) - rhs
+    auto cor = makeVecLHS();
+    m_gmres.solve(cor, res, a_tol_rel, a_tol_abs); // L(cor) = res
+    increment(a_sol, cor, RT(-1)); // sol = sol - cor
 }
 
 using GMRESPOISSON = GMRESPOISSONT<MultiFab>;
diff --git a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
index 8e90595a..da66ba00 100644
--- a/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
+++ b/ExampleCodes/LinearSolvers/GMRES/Poisson/main.cpp
@@ -25,15 +25,6 @@ int main (int argc, char* argv[])
     // size of each box (or grid)
     int max_grid_size;
 
-    // total steps in simulation
-    int nsteps;
-
-    // how often to write a plotfile
-    int plot_int;
-
-    // time step
-    amrex::Real dt;
-
     // **********************************
     // READ PARAMETER VALUES FROM INPUT DATA
     // **********************************
@@ -78,7 +69,7 @@ int main (int argc, char* argv[])
 
     // This defines the physical box, [0,1] in each direction.
     amrex::RealBox real_box({ 0., 0., 0.},
-                     { 1., 1., 1.});
+                            { 1., 1., 1.});
 
     // periodic in all direction
     amrex::Array<int,3> is_periodic{1,1,1};
@@ -121,16 +112,14 @@ int main (int argc, char* argv[])
         });
     }
 
-    // **********************************
-    // WRITE INITIAL PLOT FILE
-    // **********************************
-
-    // Write a plotfile of the initial data if plot_int > 0
     WriteSingleLevelPlotfile("rhs", rhs, {"rhs"}, geom, 0., 0);
 
-    amrex::GMRESPOISSON gmres_poisson(ba,dm);
+    amrex::GMRESPOISSON gmres_poisson(ba,dm,geom);
 
+    gmres_poisson.setVerbose(2);
     gmres_poisson.solve(phi, rhs, 1.e-12, 0.);
+
+    WriteSingleLevelPlotfile("phi", phi, {"phi"}, geom, 0., 0);
     
     }
     amrex::Finalize();