Further work on fast diagonalization

There are lots of kernels left to work on, but this is kind of the
full minimal set; ie, for one device configuration I have the full set.

fem/CMakeLists.txt

@@ -32,6 +32,7 @@ set(SRCS
   nonlininteg.cpp
   staticcond.cpp
   tmop.cpp
+  elementwisesmoother.cpp
   )
 
 set(HDRS
@@ -64,6 +65,7 @@ set(HDRS
   tfespace.hpp
   tintrules.hpp
   tmop.hpp
+  elementwisesmoother.hpp
   )
 
 if (MFEM_USE_SIDRE)

fem/bilinearform.cpp

@@ -299,6 +299,7 @@ void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
    }
 }
 
+
 void BilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
 {
    if (bbfi.Size())
@@ -326,6 +327,27 @@ void BilinearForm::AssembleElementMatrix(
    AssembleElementMatrix(i, elmat, vdofs, skip_zeros);
 }
 
+
+/// in is global, out is local to the element
+/// (this will set the size of out_element)
+void BilinearForm::ElementMatrixMult(int i, const Vector& in, Vector& out_element)
+{
+   Array<int> dofs;
+   FESpace()->GetElementDofs(i, dofs);
+   out_element.SetSize(dofs.Size());
+   out_element = 0.0;
+   if (ext)
+   {
+      ext->ElementMatrixMult(i, in, out_element);
+      return;
+   }
+   DenseMatrix elmat;
+   ComputeElementMatrix(i, elmat);
+   Vector subvec;
+   in.GetSubVector(dofs, subvec);
+   elmat.Mult(subvec, out_element);
+}
+
 void BilinearForm::AssembleElementMatrix(
    int i, const DenseMatrix &elmat, Array<int> &vdofs, int skip_zeros)
 {

fem/bilinearform.hpp

@@ -465,6 +465,14 @@ public:
    void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat,
                                  Array<int> &vdofs, int skip_zeros = 1);
 
+   /**
+      Multiply element-matrix by vector.
+
+      in is an l-vector (local to processor, but assembled across
+      elements), while out_element is local to the particular element.
+   */
+   void ElementMatrixMult(int i, const Vector& in, Vector& out_element);
+
    /// Eliminate essential boundary DOFs from the system.
    /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
        the essential part of the boundary. By default, the diagonal at the

fem/bilinearform_ext.cpp

@@ -173,4 +173,46 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
    }
 }
 
+/// TODO this needs some work
+void PABilinearFormExtension::ElementMatrixMult(int i, const Vector &x,
+                                                Vector &y_element)
+{
+   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
+   Vector yelem_preorder(y_element);
+
+   const int iSz = integrators.Size();
+   if (elem_restrict_lex)
+   {
+      elem_restrict_lex->Mult(x, localX);
+      // localY = 0.0;
+      for (int j = 0; j < iSz; ++j)
+      {
+         integrators[j]->AddMultElementPA(i, localX, yelem_preorder);
+      }
+
+      // reorder... (kinda like the transpose of elem_restrict_lex)
+      const FiniteElement *fe = trialFes->GetFE(i);
+      const TensorBasisElement* el =
+         dynamic_cast<const TensorBasisElement*>(fe);
+      MFEM_VERIFY(el, "Not a tensor element!");
+      const Array<int> &fe_dof_map = el->GetDofMap();
+      MFEM_VERIFY(fe_dof_map.Size() == yelem_preorder.Size(),
+                  "Sizes don't match!");
+      for (int j = 0; j < yelem_preorder.Size(); ++j)
+      {
+         // y_element[j] = yelem_preorder[fe_dof_map[j]];
+         y_element[fe_dof_map[j]] = yelem_preorder[j];
+      }
+
+   }
+   else
+   {
+      y_element = 0.0; // TODO: do I want this?
+      for (int i = 0; i < iSz; ++i)
+      {
+         integrators[i]->AddMultElementPA(i, x, y_element);
+      }
+   }
+}
+
 } // namespace mfem

fem/bilinearform_ext.hpp

@@ -53,6 +53,13 @@ public:
                                  OperatorHandle &A, Vector &X, Vector &B,
                                  int copy_interior = 0) = 0;
    virtual void Update() = 0;
+
+   /// Here x is local-to-processor but contains all elements
+   /// (it's an l-vector), while y_element is striclty local to an element
+   virtual void ElementMatrixMult(int i, const Vector &x, Vector &y_element)
+   {
+      mfem_error("Not implemented!");
+   }
 };
 
 /// Data and methods for fully-assembled bilinear forms
@@ -117,6 +124,7 @@ public:
    void Mult(const Vector &x, Vector &y) const;
    void MultTranspose(const Vector &x, Vector &y) const;
    void Update();
+   void ElementMatrixMult(int i, const Vector &x, Vector &y_element);
 };
 
 /// Data and methods for matrix-free bilinear forms

fem/bilininteg.cpp

@@ -34,13 +34,20 @@ void BilinearFormIntegrator::AssembleDiagonalPA(Vector&) const
 
 void BilinearFormIntegrator::AddMultPA(const Vector &, Vector &) const
 {
-   mfem_error ("BilinearFormIntegrator::MultAssembled (...)\n"
+   mfem_error ("BilinearFormIntegrator::AddMultPA (...)\n"
+               "   is not implemented for this class.");
+}
+
+void BilinearFormIntegrator::AddMultElementPA(int element,
+                                              const Vector &, Vector &) const
+{
+   mfem_error ("BilinearFormIntegrator::AddMultElementPA (...)\n"
                "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AddMultTransposePA(const Vector &, Vector &) const
 {
-   mfem_error ("BilinearFormIntegrator::MultAssembledTranspose (...)\n"
+   mfem_error ("BilinearFormIntegrator::AddMultTransposePA (...)\n"
                "   is not implemented for this class.");
 }
 

fem/bilininteg.hpp

@@ -56,6 +56,9 @@ public:
        called. */
    virtual void AddMultPA(const Vector &x, Vector &y) const;
 
+   /// Method for partially assembled action of an element matrix.
+   virtual void AddMultElementPA(int element, const Vector &x, Vector &y) const;
+
    /// Method for partially assembled transposed action.
    /** Perform the transpose action of integrator on the input @a x and add the
        result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
@@ -1733,6 +1736,8 @@ public:
 
    virtual void AddMultPA(const Vector&, Vector&) const;
 
+   virtual void AddMultElementPA(int element, const Vector &x, Vector &y) const;
+
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe);
 
@@ -1777,6 +1782,8 @@ public:
 
    virtual void AddMultPA(const Vector&, Vector&) const;
 
+   virtual void AddMultElementPA(int element, const Vector &x, Vector &y) const;
+   
    virtual MassIntegrator* Copy() const;
 
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,

fem/bilininteg_mass.cpp

@@ -121,7 +121,6 @@ static void PAMassAssembleDiagonal2D(const int NE,
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
    auto B = Reshape(_B.Read(), Q1D, D1D);
-   // auto Bt = Reshape(_Bt.Read(), D1D, Q1D); // ?? (TODO atb@llnl.gov)
    auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
    auto y = Reshape(_diag.ReadWrite(), D1D, D1D, NE);
    MFEM_FORALL(e, NE,
@@ -173,7 +172,6 @@ static void PAMassAssembleDiagonal3D(const int NE,
    MFEM_VERIFY(D1D <= MAX_D1D, "");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "");
    auto B = Reshape(_B.Read(), Q1D, D1D);
-   // auto Bt = Reshape(_Bt.Read(), D1D, Q1D); // ?? (TODO atb@llnl.gov)
    auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
    auto y = Reshape(_diag.ReadWrite(), D1D, D1D, D1D, NE);
    MFEM_FORALL(e, NE,
@@ -470,6 +468,117 @@ static void PAMassApply2D(const int NE,
    });
 }
 
+template<const int T_D1D = 0,
+         const int T_Q1D = 0,
+         const int T_NBZ = 0>
+static void SmemPAMassApply2DElement(int e,
+                                     DeviceTensor<2,const double> b,
+                                     DeviceTensor<3,const double> op,
+                                     DeviceTensor<2,const double> x,
+                                     DeviceTensor<2,double> y,
+                                     const int d1d = 0,
+                                     const int q1d = 0)
+{
+   const int tidz = MFEM_THREAD_ID(z);
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   constexpr int NBZ = T_NBZ ? T_NBZ : 1;
+   constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
+   constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
+   constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
+   MFEM_SHARED double BBt[MQ1*MD1];
+   double (*B)[MD1] = (double (*)[MD1]) BBt;
+   double (*Bt)[MQ1] = (double (*)[MQ1]) BBt;
+   MFEM_SHARED double sm0[NBZ][MDQ*MDQ];
+   MFEM_SHARED double sm1[NBZ][MDQ*MDQ];
+   double (*X)[MD1] = (double (*)[MD1]) (sm0 + tidz);
+   double (*DQ)[MQ1] = (double (*)[MQ1]) (sm1 + tidz);
+   double (*QQ)[MQ1] = (double (*)[MQ1]) (sm0 + tidz);
+   double (*QD)[MD1] = (double (*)[MD1]) (sm1 + tidz);
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         // X[dy][dx] = x(dx,dy,e);
+         X[dy][dx] = x(dx,dy);
+      }
+   }
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(d,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(q,x,Q1D)
+         {
+            B[q][d] = b(q,d);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double dq = 0.0;
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            dq += X[dy][dx] * B[qx][dx];
+         }
+         DQ[dy][qx] = dq;
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double qq = 0.0;
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            qq += DQ[dy][qx] * B[qy][dy];
+         }
+         QQ[qy][qx] = qq * op(qx, qy, e);
+      }
+   }
+   MFEM_SYNC_THREAD;
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(d,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(q,x,Q1D)
+         {
+            Bt[d][q] = b(q,d);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double dq = 0.0;
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            dq += QQ[qy][qx] * Bt[dx][qx];
+         }
+         QD[qy][dx] = dq;
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double dd = 0.0;
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            dd += (QD[qy][dx] * Bt[dy][qy]);
+         }
+         // y(dx, dy, e) += dd;
+         y(dx, dy) += dd;
+      }
+   }
+}
+
 template<const int T_D1D = 0,
          const int T_Q1D = 0,
          const int T_NBZ = 0>
@@ -495,102 +604,15 @@ static void SmemPAMassApply2D(const int NE,
    auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
    MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
    {
-      const int tidz = MFEM_THREAD_ID(z);
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int NBZ = T_NBZ ? T_NBZ : 1;
-      constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
-      constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
-      constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
-      MFEM_SHARED double BBt[MQ1*MD1];
-      double (*B)[MD1] = (double (*)[MD1]) BBt;
-      double (*Bt)[MQ1] = (double (*)[MQ1]) BBt;
-      MFEM_SHARED double sm0[NBZ][MDQ*MDQ];
-      MFEM_SHARED double sm1[NBZ][MDQ*MDQ];
-      double (*X)[MD1] = (double (*)[MD1]) (sm0 + tidz);
-      double (*DQ)[MQ1] = (double (*)[MQ1]) (sm1 + tidz);
-      double (*QQ)[MQ1] = (double (*)[MQ1]) (sm0 + tidz);
-      double (*QD)[MD1] = (double (*)[MD1]) (sm1 + tidz);
-      MFEM_FOREACH_THREAD(dy,y,D1D)
-      {
-         MFEM_FOREACH_THREAD(dx,x,D1D)
-         {
-            X[dy][dx] = x(dx,dy,e);
-         }
-      }
-      if (tidz == 0)
-      {
-         MFEM_FOREACH_THREAD(d,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(q,x,Q1D)
-            {
-               B[q][d] = b(q,d);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(dy,y,D1D)
-      {
-         MFEM_FOREACH_THREAD(qx,x,Q1D)
-         {
-            double dq = 0.0;
-            for (int dx = 0; dx < D1D; ++dx)
-            {
-               dq += X[dy][dx] * B[qx][dx];
-            }
-            DQ[dy][qx] = dq;
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(qy,y,Q1D)
-      {
-         MFEM_FOREACH_THREAD(qx,x,Q1D)
-         {
-            double qq = 0.0;
-            for (int dy = 0; dy < D1D; ++dy)
-            {
-               qq += DQ[dy][qx] * B[qy][dy];
-            }
-            QQ[qy][qx] = qq * op(qx, qy, e);
-         }
-      }
-      MFEM_SYNC_THREAD;
-      if (tidz == 0)
-      {
-         MFEM_FOREACH_THREAD(d,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(q,x,Q1D)
-            {
-               Bt[d][q] = b(q,d);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(qy,y,Q1D)
-      {
-         MFEM_FOREACH_THREAD(dx,x,D1D)
-         {
-            double dq = 0.0;
-            for (int qx = 0; qx < Q1D; ++qx)
-            {
-               dq += QQ[qy][qx] * Bt[dx][qx];
-            }
-            QD[qy][dx] = dq;
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(dy,y,D1D)
-      {
-         MFEM_FOREACH_THREAD(dx,x,D1D)
-         {
-            double dd = 0.0;
-            for (int qy = 0; qy < Q1D; ++qy)
-            {
-               dd += (QD[qy][dx] * Bt[dy][qy]);
-            }
-            y(dx, dy, e) += dd;
-         }
-      }
+      auto x_element = DeviceTensor<2, const double>(
+         (const double*) x + e * D1D * D1D,
+         D1D, D1D);
+      auto y_element = DeviceTensor<2, double>(
+         (double*) y + e * D1D * D1D,
+         D1D, D1D);
+
+      SmemPAMassApply2DElement<T_D1D, T_Q1D, T_NBZ>(e, b, op, x_element,
+                                                    y_element, d1d, q1d);
    });
 }
 
@@ -736,14 +758,174 @@ static void PAMassApply3D(const int NE,
    });
 }
 
+
+template<const int T_D1D = 0,
+         const int T_Q1D = 0>
+static void SmemPAMassApply3DElement(int e,
+                                     DeviceTensor<2,const double> b,
+                                     DeviceTensor<4,const double> op,
+                                     DeviceTensor<3,const double> x,
+                                     DeviceTensor<3,double> y,
+                                     const int d1d = 0,
+                                     const int q1d = 0)
+{
+   const int tidz = MFEM_THREAD_ID(z);
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
+   constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
+   constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
+   MFEM_SHARED double sDQ[MQ1*MD1];
+   double (*B)[MD1] = (double (*)[MD1]) sDQ;
+   double (*Bt)[MQ1] = (double (*)[MQ1]) sDQ;
+   MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
+   MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
+   double (*X)[MD1][MD1]   = (double (*)[MD1][MD1]) sm0;
+   double (*DDQ)[MD1][MQ1] = (double (*)[MD1][MQ1]) sm1;
+   double (*DQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm0;
+   double (*QQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm1;
+   double (*QQD)[MQ1][MD1] = (double (*)[MQ1][MD1]) sm0;
+   double (*QDD)[MD1][MD1] = (double (*)[MD1][MD1]) sm1;
+   MFEM_FOREACH_THREAD(dz,z,D1D)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(dx,x,D1D)
+         {
+            // X[dz][dy][dx] = x(dx,dy,dz,e);
+            X[dz][dy][dx] = x(dx,dy,dz);
+         }
+      }
+   }
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(d,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(q,x,Q1D)
+         {
+            B[q][d] = b(q,d);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dz,z,D1D)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            double u = 0.0;
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               u += X[dz][dy][dx] * B[qx][dx];
+            }
+            DDQ[dz][dy][qx] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dz,z,D1D)
+   {
+      MFEM_FOREACH_THREAD(qy,y,Q1D)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            double u = 0.0;
+            for (int dy = 0; dy < D1D; ++dy)
+            {
+               u += DDQ[dz][dy][qx] * B[qy][dy];
+            }
+            DQQ[dz][qy][qx] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qz,z,Q1D)
+   {
+      MFEM_FOREACH_THREAD(qy,y,Q1D)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            double u = 0.0;
+            for (int dz = 0; dz < D1D; ++dz)
+            {
+               u += DQQ[dz][qy][qx] * B[qz][dz];
+            }
+            QQQ[qz][qy][qx] = u * op(qx,qy,qz,e);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(d,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(q,x,Q1D)
+         {
+            Bt[d][q] = b(q,d);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qz,z,Q1D)
+   {
+      MFEM_FOREACH_THREAD(qy,y,Q1D)
+      {
+         MFEM_FOREACH_THREAD(dx,x,D1D)
+         {
+            double u = 0.0;
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               u += QQQ[qz][qy][qx] * Bt[dx][qx];
+            }
+            QQD[qz][qy][dx] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qz,z,Q1D)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(dx,x,D1D)
+         {
+            double u = 0.0;
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               u += QQD[qz][qy][dx] * Bt[dy][qy];
+            }
+            QDD[qz][dy][dx] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dz,z,D1D)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(dx,x,D1D)
+         {
+            double u = 0.0;
+            for (int qz = 0; qz < Q1D; ++qz)
+            {
+               u += QDD[qz][dy][dx] * Bt[dz][qz];
+            }
+            // y(dx,dy,dz,e) += u;
+            y(dx,dy,dz) += u;
+         }
+      }
+   }
+}
+
+
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
 static void SmemPAMassApply3D(const int NE,
-                              const Array<double> &b_,
-                              const Array<double> &bt_,
-                              const Vector &op_,
-                              const Vector &x_,
-                              Vector &y_,
+                              const Array<double> &_b,
+                              const Array<double> &_bt,
+                              const Vector &_op,
+                              const Vector &_x,
+                              Vector &_y,
                               const int d1d = 0,
                               const int q1d = 0)
 {
@@ -753,156 +935,22 @@ static void SmemPAMassApply3D(const int NE,
    constexpr int M1D = T_D1D ? T_D1D : MAX_D1D;
    MFEM_VERIFY(D1D <= M1D, "");
    MFEM_VERIFY(Q1D <= M1Q, "");
-   auto b = Reshape(b_.Read(), Q1D, D1D);
-   auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
-   auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
-   auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
+   auto b = Reshape(_b.Read(), Q1D, D1D);
+   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
+   auto x = Reshape(_x.Read(), D1D, D1D, D1D, NE);
+   auto y = Reshape(_y.ReadWrite(), D1D, D1D, D1D, NE);
    MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
    {
-      const int tidz = MFEM_THREAD_ID(z);
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
-      constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
-      constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
-      MFEM_SHARED double sDQ[MQ1*MD1];
-      double (*B)[MD1] = (double (*)[MD1]) sDQ;
-      double (*Bt)[MQ1] = (double (*)[MQ1]) sDQ;
-      MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
-      MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
-      double (*X)[MD1][MD1]   = (double (*)[MD1][MD1]) sm0;
-      double (*DDQ)[MD1][MQ1] = (double (*)[MD1][MQ1]) sm1;
-      double (*DQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm0;
-      double (*QQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm1;
-      double (*QQD)[MQ1][MD1] = (double (*)[MQ1][MD1]) sm0;
-      double (*QDD)[MD1][MD1] = (double (*)[MD1][MD1]) sm1;
-      MFEM_FOREACH_THREAD(dz,z,D1D)
-      {
-         MFEM_FOREACH_THREAD(dy,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(dx,x,D1D)
-            {
-               X[dz][dy][dx] = x(dx,dy,dz,e);
-            }
-         }
-      }
-      if (tidz == 0)
-      {
-         MFEM_FOREACH_THREAD(d,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(q,x,Q1D)
-            {
-               B[q][d] = b(q,d);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(dz,z,D1D)
-      {
-         MFEM_FOREACH_THREAD(dy,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(qx,x,Q1D)
-            {
-               double u = 0.0;
-               for (int dx = 0; dx < D1D; ++dx)
-               {
-                  u += X[dz][dy][dx] * B[qx][dx];
-               }
-               DDQ[dz][dy][qx] = u;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(dz,z,D1D)
-      {
-         MFEM_FOREACH_THREAD(qy,y,Q1D)
-         {
-            MFEM_FOREACH_THREAD(qx,x,Q1D)
-            {
-               double u = 0.0;
-               for (int dy = 0; dy < D1D; ++dy)
-               {
-                  u += DDQ[dz][dy][qx] * B[qy][dy];
-               }
-               DQQ[dz][qy][qx] = u;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(qz,z,Q1D)
-      {
-         MFEM_FOREACH_THREAD(qy,y,Q1D)
-         {
-            MFEM_FOREACH_THREAD(qx,x,Q1D)
-            {
-               double u = 0.0;
-               for (int dz = 0; dz < D1D; ++dz)
-               {
-                  u += DQQ[dz][qy][qx] * B[qz][dz];
-               }
-               QQQ[qz][qy][qx] = u * op(qx,qy,qz,e);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      if (tidz == 0)
-      {
-         MFEM_FOREACH_THREAD(d,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(q,x,Q1D)
-            {
-               Bt[d][q] = b(q,d);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(qz,z,Q1D)
-      {
-         MFEM_FOREACH_THREAD(qy,y,Q1D)
-         {
-            MFEM_FOREACH_THREAD(dx,x,D1D)
-            {
-               double u = 0.0;
-               for (int qx = 0; qx < Q1D; ++qx)
-               {
-                  u += QQQ[qz][qy][qx] * Bt[dx][qx];
-               }
-               QQD[qz][qy][dx] = u;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(qz,z,Q1D)
-      {
-         MFEM_FOREACH_THREAD(dy,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(dx,x,D1D)
-            {
-               double u = 0.0;
-               for (int qy = 0; qy < Q1D; ++qy)
-               {
-                  u += QQD[qz][qy][dx] * Bt[dy][qy];
-               }
-               QDD[qz][dy][dx] = u;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(dz,z,D1D)
-      {
-         MFEM_FOREACH_THREAD(dy,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(dx,x,D1D)
-            {
-               double u = 0.0;
-               for (int qz = 0; qz < Q1D; ++qz)
-               {
-                  u += QDD[qz][dy][dx] * Bt[dz][qz];
-               }
-               y(dx,dy,dz,e) += u;
-            }
-         }
-      }
+      auto x_element = DeviceTensor<3, const double>(
+         (const double*) x + e * D1D * D1D * D1D,
+         D1D, D1D, D1D);
+      auto y_element = DeviceTensor<3, double>(
+         (double*) y + e * D1D * D1D * D1D,
+         D1D, D1D, D1D);
+
+      SmemPAMassApply3DElement<T_D1D, T_Q1D>(e, b, op, x_element,
+                                             y_element, D1D, Q1D);
+
    });
 }
 
@@ -969,6 +1017,61 @@ void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    PAMassApply(dim, dofs1D, quad1D, ne, maps->B, maps->Bt, pa_data, x, y);
 }
 
+/**
+   as written this wants a global e-vector for x and a local
+   element vector for y
+
+   not sure that's the right interface
+*/
+void MassIntegrator::AddMultElementPA(int element, const Vector &x,
+                                      Vector &y) const
+{
+   if (dim == 2)
+   {
+      auto B = Reshape(maps->B.Read(), quad1D, dofs1D);
+      auto op = Reshape(pa_data.Read(), quad1D, quad1D, ne);
+
+      auto _x = Reshape(x.Read(), dofs1D, dofs1D, ne);
+      auto _y = Reshape(y.ReadWrite(), dofs1D, dofs1D);
+
+      auto x_element = DeviceTensor<2, const double>(
+         (const double*) _x + element * dofs1D * dofs1D,
+         dofs1D, dofs1D);
+
+      SmemPAMassApply2DElement(element,
+                               B,
+                               op,
+                               x_element,
+                               _y,
+                               dofs1D,
+                               quad1D);
+   }
+   else if (dim == 3)
+   {
+      auto B = Reshape(maps->B.Read(), quad1D, dofs1D);
+      auto op = Reshape(pa_data.Read(), quad1D, quad1D, quad1D, ne);
+      auto _x = Reshape(x.Read(), dofs1D, dofs1D, dofs1D, ne);
+      auto _y = Reshape(y.ReadWrite(), dofs1D, dofs1D, dofs1D);
+
+      auto x_element = DeviceTensor<3, const double>(
+         (const double*) _x + element * dofs1D * dofs1D * dofs1D,
+         dofs1D, dofs1D, dofs1D);
+
+      SmemPAMassApply3DElement(element,
+                               B,
+                               op,
+                               x_element,
+                               _y,
+                               dofs1D,
+                               quad1D);
+   }
+   else
+   {
+      mfem_error("Not implemented!");
+   }
+
+}
+
 MassIntegrator* MassIntegrator::Copy() const
 {
    return new MassIntegrator(*this);

fem/elementwisesmoother.cpp

@@ -0,0 +1,526 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "fem.hpp"
+#include "../general/forall.hpp"
+
+namespace mfem
+{
+
+ElementWiseSmoother::ElementWiseSmoother(const mfem::FiniteElementSpace& fespace)
+   :
+   mfem::Solver(fespace.GetVSize()),
+   fespace_(fespace)
+{
+}
+
+void ElementWiseSmoother::Mult(const Vector& b, Vector& x) const
+{
+   Vector local_residual;
+   Vector local_correction;
+   Vector b_local;
+   Array<int> local_dofs;
+
+   if (!iterative_mode)
+   {
+      x = 0.0;
+   }
+
+   for (int e = 0; e < fespace_.GetNE(); ++e)
+   {
+      fespace_.GetElementDofs(e, local_dofs);
+      b.GetSubVector(local_dofs, b_local);
+      local_correction.SetSize(b_local.Size());
+      ElementResidual(e, b_local, x, local_residual);
+
+      LocalSmoother(e, local_residual, local_correction);
+
+      x.AddElementVector(local_dofs, local_correction);
+   }
+}
+
+void ElementWiseSmoother::ElementResidual(int e, const Vector& b_local,
+                                          const Vector& x,
+                                          Vector& r_local) const
+{
+   r_local.SetSize(b_local.Size());
+   r_local = 0.0;
+
+   GetElementFromMatVec(e, x, r_local);
+   r_local -= b_local;
+   r_local *= -1.0;
+}
+
+AdditiveSchwarzLORSmoother::AdditiveSchwarzLORSmoother(const mfem::FiniteElementSpace& fespace,
+                                                   const Array<int>& ess_tdof_list,
+                                                   mfem::BilinearForm& aform,
+                                                   const mfem::Vector& diag,
+                                                   mfem::SparseMatrix* LORmat,
+                                                   double scale)
+   :
+   mfem::Solver(fespace.GetVSize()),
+   fespace_(fespace),
+   ess_tdof_list_(ess_tdof_list),
+   aform_(aform),
+   diag_(diag),
+   LORmat_(LORmat),
+   scale_(scale)
+{
+   const Table& el_dof = fespace_.GetElementToDofTable();
+   Table dof_el;
+   mfem::Transpose(el_dof, dof_el);
+   mfem::Mult(el_dof, dof_el, el_to_el_);
+
+   countingVector.SetSize(fespace_.GetTrueVSize());
+   countingVector = 0.0;
+
+   Array<int> local_dofs;
+
+   for (int e = 0; e < fespace_.GetNE(); ++e)
+   {
+      fespace_.GetElementDofs(e, local_dofs);
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         countingVector[local_dofs[i]] += 1.0;
+      }
+   }
+
+   for (int i = 0; i < fespace_.GetTrueVSize(); ++i)
+   {
+      countingVector[i] = 1.0 / std::sqrt(countingVector[i]);
+   }
+}
+
+void AdditiveSchwarzLORSmoother::Mult(const Vector& b, Vector& x) const
+{
+   x = 0.0;
+   Array<int> local_dofs;
+   Vector b_local;
+   Vector x_local;
+
+   for (int e = 0; e < fespace_.GetNE(); ++e)
+   {
+      fespace_.GetElementDofs(e, local_dofs);
+      b.GetSubVector(local_dofs, b_local);
+      x_local.SetSize(b_local.Size());
+
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         b_local[i] *= countingVector[local_dofs[i]];
+      }
+
+      LocalSmoother(e, b_local, x_local);
+
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         x_local[i] *= countingVector[local_dofs[i]];
+      }
+
+      x.AddElementVector(local_dofs, x_local);
+   }
+
+   x *= scale_;
+
+   auto I = ess_tdof_list_.Read();
+   auto B = b.Read();
+   auto X = x.Write();
+   MFEM_FORALL(i, ess_tdof_list_.Size(), X[I[i]] = B[I[i]]; );
+}
+
+void AdditiveSchwarzLORSmoother::LocalSmoother(int e, const Vector& in, Vector& out) const
+{
+   Array<int> local_dofs;
+   fespace_.GetElementDofs(e, local_dofs);
+   DenseMatrix elmat(local_dofs.Size(), local_dofs.Size());
+   LORmat_->GetSubMatrix(local_dofs, local_dofs, elmat);
+   DenseMatrixInverse inv(elmat);
+   inv.Mult(in, out);
+}
+
+AdditiveSchwarzApproxLORSmoother::AdditiveSchwarzApproxLORSmoother(const mfem::FiniteElementSpace& fespace,
+                                                   const Array<int>& ess_tdof_list,
+                                                   mfem::BilinearForm& aform,
+                                                   const mfem::Vector& diag,
+                                                   const mfem::Vector& LORdiag,
+                                                   mfem::SparseMatrix* LORmat,
+                                                   double scale)
+   :
+   mfem::Solver(fespace.GetVSize()),
+   fespace_(fespace),
+   ess_tdof_list_(ess_tdof_list),
+   aform_(aform),
+   diag_(diag),
+   LORdiag_(LORdiag),
+   LORmat_(LORmat),
+   scale_(scale)
+{
+   const Table& el_dof = fespace_.GetElementToDofTable();
+   Table dof_el;
+   mfem::Transpose(el_dof, dof_el);
+   mfem::Mult(el_dof, dof_el, el_to_el_);
+
+   countingVector.SetSize(fespace_.GetTrueVSize());
+   countingVector = 0.0;
+
+   Array<int> local_dofs;
+   Array<int> local_dofs_lex;
+
+   for (int e = 0; e < fespace_.GetNE(); ++e)
+   {
+      fespace_.GetElementDofs(e, local_dofs);
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         countingVector[local_dofs[i]] += 1.0;
+      }
+
+      Array<int> neighbors;
+      el_to_el_.GetRow(e, neighbors);
+   }
+
+   for (int i = 0; i < fespace_.GetTrueVSize(); ++i)
+   {
+      countingVector[i] = 1.0 / std::sqrt(countingVector[i]);
+   }
+
+   A1.resize(fespace_.GetNE());
+   A2.resize(fespace_.GetNE());
+   B1.resize(fespace_.GetNE());
+   B2.resize(fespace_.GetNE());
+   invA2.resize(fespace_.GetNE());
+   invB1.resize(fespace_.GetNE());
+   schurA1.resize(fespace_.GetNE());
+   schurB2.resize(fespace_.GetNE());
+   syl.resize(fespace_.GetNE());
+   inv.resize(fespace_.GetNE());
+
+   for (int e = 0; e < fespace_.GetNE(); ++e)
+   {
+      const FiniteElement& el = *fespace_.GetFE(e);
+      const TensorBasisElement* ltel =
+         dynamic_cast<const TensorBasisElement*>(&el);
+      MFEM_VERIFY(ltel, "FE space must be tensor product space");
+      lexdofmap = ltel->GetDofMap();
+
+      fespace_.GetElementDofs(e, local_dofs);
+      elmat = DenseMatrix(local_dofs.Size(), local_dofs.Size());
+
+      local_dofs_lex.SetSize(local_dofs.Size());
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         local_dofs_lex[i] = local_dofs[lexdofmap[i]];
+      }
+
+      LORmat_->GetSubMatrix(local_dofs_lex, local_dofs_lex, elmat);
+      inv[e].SetOperator(elmat);
+
+      int N1D = std::sqrt(local_dofs_lex.Size());
+      DenseMatrix Atilde = elmat;
+
+      for (int i = 0; i < N1D; ++i)
+      {
+         for (int j = 0; j < N1D; ++j)
+         {
+            DenseMatrix blockMat(N1D, N1D);
+            for (int k = 0; k < N1D; ++k)
+            {
+               for (int l = 0; l < N1D; ++l)
+               {
+                  blockMat(k,l) = elmat(k + i*N1D, l + j*N1D);
+               }
+            }
+
+            Atilde.SetRow(j + i * N1D, Vector(blockMat.Data(), N1D * N1D));
+         }
+      }
+
+      A1[e].SetSize(N1D, N1D);
+      A2[e].SetSize(N1D, N1D);
+      B1[e].SetSize(N1D, N1D);
+      B2[e].SetSize(N1D, N1D);
+
+      DenseMatrixSVD svd(Atilde);
+      svd.Eval(Atilde);
+
+      // for (int i = 0; i < local_dofs.Size(); ++i)
+      // {
+      //    std::cout << svd.Singularvalue(i) << " ";
+      // }
+      // std::cout << std::endl;
+
+      Vector A1vec(A1[e].Data(), N1D * N1D);
+      svd.GetU().GetColumn(0, A1vec);
+      A1[e] *= std::sqrt(svd.Singularvalue(0));
+
+      Vector A2vec(A2[e].Data(), N1D * N1D);
+      svd.GetU().GetColumn(1, A2vec);
+      A2[e] *= std::sqrt(svd.Singularvalue(1));
+
+      Vector B1vec(B1[e].Data(), N1D * N1D);
+      svd.GetV_T().GetRow(0, B1vec);
+      B1[e] *= std::sqrt(svd.Singularvalue(0));
+
+      Vector B2vec(B2[e].Data(), N1D * N1D);
+      svd.GetV_T().GetRow(1, B2vec);
+      B2[e] *= std::sqrt(svd.Singularvalue(1));
+
+
+      ////---- Test if "tensorization" worked
+
+      // Vector xTest(local_dofs.Size());
+      // Vector yTest(local_dofs.Size());
+      // Vector yTest2(local_dofs.Size());
+      // Vector yTest3(local_dofs.Size());
+
+      // xTest.Randomize();
+
+      // tic_toc.Clear();
+      // tic_toc.Start();
+      // elmat.Mult(xTest, yTest);
+      // tic_toc.Stop();
+      // std::cout << "elmat mult: " << tic_toc.RealTime() << std::endl;
+
+      // tic_toc.Clear();
+      // tic_toc.Start();
+      // TensorProductMult2D(A1[e], B1[e], xTest, yTest2);
+      // TensorProductMult2D(A2[e], B2[e], xTest, yTest3);
+      // yTest2 += yTest3;
+      // tic_toc.Stop();
+      // std::cout << "tp mult:    " << tic_toc.RealTime() << std::endl;
+
+      // yTest -= yTest2;
+
+      // std::cout << "norm = " << yTest.Norml2() << std::endl;
+
+      ////---- Test end
+
+      invA2[e].SetOperator(A2[e]);
+      invB1[e].SetOperator(B1[e]);
+
+      invA2[e].Mult(A1[e]); // A1 = inv(A2) * A1
+      invB1[e].Mult(B2[e]); // B2 = inv(B1) * B2
+
+      schurA1[e] = new DenseMatrixSchurDecomposition(A1[e]);
+      schurB2[e] = new DenseMatrixSchurDecomposition(B2[e]);
+
+      // ////---- Test if Schur decomposition worked
+      // DenseMatrix testMat(A1.Height(), A1.Width());
+      // mfem::Mult(schurA1->GetT(), schurA1->GetQ_T(), testMat);
+      // DenseMatrix tmp(testMat);
+      // mfem::Mult(schurA1->GetQ(), tmp, testMat);
+
+      // std::cout << "A1 = " << A1.MaxMaxNorm() << std::endl;
+      // std::cout << "testMat = " << testMat.MaxMaxNorm() << std::endl;
+
+      // testMat -= A1;
+      // std::cout << "diff = " << testMat.MaxMaxNorm() << std::endl;
+      // ////---- Test end
+
+      syl[e] = new DenseMatrixSylvesterSolver(schurB2[e]->GetT(), schurA1[e]->GetT(), false, true);
+   }
+}
+
+void AdditiveSchwarzApproxLORSmoother::Mult(const Vector& b, Vector& x) const
+{
+   x = 0.0;
+   Array<int> local_dofs;
+   Vector b_local, b_local_lex;
+   Vector x_local, x_local_lex;
+
+   for (int e = 0; e < fespace_.GetNE(); ++e)
+   {
+      fespace_.GetElementDofs(e, local_dofs);
+      b.GetSubVector(local_dofs, b_local);
+      x_local.SetSize(b_local.Size());
+      b_local_lex.SetSize(b_local.Size());
+      x_local_lex.SetSize(b_local.Size());
+
+      double minval = 1e300;
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         minval = std::min(minval, diag_[local_dofs[i]]);
+      }
+
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         b_local[i] *= countingVector[local_dofs[i]];
+         // b_local[i] *= diag_[local_dofs[i]];
+         // b_local[i] *= minval;
+         // b_local[i] *= sqrt(elmat(i,i) / LORdiag_[local_dofs[i]]);
+      }
+
+      for (int i = 0; i < lexdofmap.Size(); ++i)
+      {
+         b_local_lex[i] = b_local[lexdofmap[i]];
+      }
+
+      LocalSmoother(e, b_local_lex, x_local_lex);
+
+      for (int i = 0; i < lexdofmap.Size(); ++i)
+      {
+         x_local[lexdofmap[i]] = x_local_lex[i];
+      }
+
+      for (int i = 0; i < local_dofs.Size(); ++i)
+      {
+         // x_local[i] *= sqrt(elmat(i,i) / LORdiag_[local_dofs[i]]);
+         // x_local[i] *= minval;
+         // x_local[i] *= diag_[local_dofs[i]];
+         x_local[i] *= countingVector[local_dofs[i]];
+      }
+
+      x.AddElementVector(local_dofs, x_local);
+   }
+
+   x *= scale_;
+
+   auto I = ess_tdof_list_.Read();
+   auto B = b.Read();
+   auto X = x.Write();
+   MFEM_FORALL(i, ess_tdof_list_.Size(), X[I[i]] = B[I[i]]; );
+}
+
+void AdditiveSchwarzApproxLORSmoother::LocalSmoother(int e, const Vector& in, Vector& out) const
+{
+   // TensorProductMult2D(invA2[e], invB1[e], in, out);
+   // TensorProductMult2D(schurA1[e]->GetQ_T(), schurB2[e]->GetQ_T(), out, out);
+   // syl[e]->Mult(out, out);
+   // TensorProductMult2D(schurA1[e]->GetQ(), schurB2[e]->GetQ(), out, out);
+   
+   inv[e].Mult(in, out);
+}
+
+void AdditiveSchwarzApproxLORSmoother::TensorProductMult2D(const DenseMatrix& A, const DenseMatrix& B, const Vector& in, Vector& out) const
+{
+   const int N1D = std::sqrt(in.Size());
+   DenseMatrix IN(in.GetData(), N1D, N1D);
+   DenseMatrix OUT(out.GetData(), N1D, N1D);
+   DenseMatrix tmp(N1D, N1D);
+
+   mfem::MultABt(IN, A, tmp);
+   mfem::Mult(B, tmp, OUT);
+}
+
+// Computes out = (A (x) B) * in
+void AdditiveSchwarzApproxLORSmoother::TensorProductMult2D(const Operator& A, const Operator& B, const Vector& in, Vector& out) const
+{
+   const int N1D = std::sqrt(in.Size());
+   Vector x1D(N1D);
+   Vector y1D(N1D);
+
+   auto in_ = Reshape(in.Read(), N1D, N1D);
+   auto out_ = Reshape(out.ReadWrite(), N1D, N1D);
+
+   for (int i = 0; i < N1D; ++i)
+   {
+      for (int j = 0; j < N1D; ++j)
+      {
+         x1D(j) = in_(j,i);
+      }
+      
+      B.Mult(x1D, y1D);
+
+      for (int j = 0; j < N1D; ++j)
+      {
+         out_(j,i) = y1D(j);
+      }
+   }
+
+   for (int i = 0; i < N1D; ++i)
+   {
+      for (int j = 0; j < N1D; ++j)
+      {
+         x1D(j) = out_(i,j);
+      }
+      
+      A.Mult(x1D, y1D);
+
+      for (int j = 0; j < N1D; ++j)
+      {
+         out_(i,j) = y1D(j);
+      }
+   }
+}
+
+
+
+
+
+ElementWiseJacobi::ElementWiseJacobi(const mfem::FiniteElementSpace& fespace,
+                                     mfem::BilinearForm& aform,
+                                     const mfem::Vector& global_diag,
+                                     double scale)
+   :
+   ElementWiseSmoother(fespace),
+   aform_(aform),
+   global_diag_(global_diag),
+   scale_(scale)
+{
+   const Table& el_dof = fespace_.GetElementToDofTable();
+   Table dof_el;
+   mfem::Transpose(el_dof, dof_el);
+   mfem::Mult(el_dof, dof_el, el_to_el_);
+}
+
+// x is *global*, y is local to the element
+void ElementWiseJacobi::GetElementFromMatVec(int element, const mfem::Vector& x,
+                                             mfem::Vector& y_element) const
+{
+   Array<int> neighbors;
+   el_to_el_.GetRow(element, neighbors);
+   Array<int> row_dofs;
+   Array<int> col_dofs;
+   fespace_.GetElementDofs(element, row_dofs);
+   Vector x_local;
+   Vector z_element;
+   DenseMatrix elmat;
+
+   y_element.SetSize(row_dofs.Size());
+   z_element.SetSize(row_dofs.Size());
+   y_element = 0.0;
+
+   for (int i = 0; i < neighbors.Size(); ++i)
+   {
+      int ne = neighbors[i];
+      fespace_.GetElementDofs(ne, col_dofs);
+
+      z_element = 0.0;
+      aform_.ElementMatrixMult(ne, x, z_element);
+
+      // okay, this next section seems very inefficient
+      // (could probably store and precompute some kind of map?)
+      for (int j = 0; j < row_dofs.Size(); ++j)
+      {
+         const int rd = row_dofs[j];
+         for (int k = 0; k < col_dofs.Size(); ++k)
+         {
+            const int cd = col_dofs[k];
+            if (rd == cd)
+            {
+               y_element[j] += z_element[k];
+               break;
+            }
+         }
+      }
+   }
+}
+
+void ElementWiseJacobi::LocalSmoother(int e, const Vector& in, Vector& out) const
+{
+   DenseMatrix elmat;
+   Array<int> local_dofs;
+   fespace_.GetElementDofs(e, local_dofs);
+   for (int i = 0; i < in.Size(); ++i)
+   {
+      out[i] = (scale_ / global_diag_(local_dofs[i])) * in[i];
+   }
+}
+
+}

fem/elementwisesmoother.hpp

@@ -0,0 +1,176 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_ELEMENTWISESMOOTHER
+#define MFEM_ELEMENTWISESMOOTHER
+
+#include "../config/config.hpp"
+#include "fespace.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+namespace mfem
+{
+
+/**
+   Applies a smoother element-by-element.
+
+   The interface basically requires you to implement GetElementFromMatVec()
+   and LocalSmoother()
+*/
+class ElementWiseSmoother : public mfem::Solver
+{
+public:
+   ElementWiseSmoother(const mfem::FiniteElementSpace& fespace);
+   virtual ~ElementWiseSmoother() {}
+
+   virtual void SetOperator(const mfem::Operator &op) {}
+
+   virtual void Mult(const mfem::Vector& in, mfem::Vector& out) const;
+
+   /**
+      b and r are local to the element, but x has to be global.
+   */
+   virtual void ElementResidual(int element, const mfem::Vector& b,
+                                const mfem::Vector& x,
+                                mfem::Vector& r) const;
+
+   /**
+      This produces the same result in y_element as if you did a
+      global matvec with (global) x, and then extracted the
+      dofs corresponding to a particular element.
+
+      Since elements have a bounded number of neighbors, this can
+      be done more efficiently than a global matvec.
+   */
+   virtual void GetElementFromMatVec(int element, const mfem::Vector& x,
+                                     mfem::Vector& y_element) const = 0;
+
+   /**
+      Interface is that `in` should be a residual, and then `out` is a correction
+      to the solution. (Ie, in is not a right-hand side, and out is not a current
+      iterate; out gets totally overwritten, not updated.)
+
+      I am not sure this is the *correct* interface for this use, but it is the
+      *current* interface.
+   */
+   virtual void LocalSmoother(int element, const mfem::Vector& in,
+                              mfem::Vector& out) const = 0;
+
+protected:
+   const mfem::FiniteElementSpace& fespace_;
+};
+
+class AdditiveSchwarzLORSmoother : public Solver
+{
+ public:
+   AdditiveSchwarzLORSmoother(const mfem::FiniteElementSpace& fespace,
+                              const Array<int>& ess_tdof_list,
+                              mfem::BilinearForm& aform,
+                              const mfem::Vector& diag,
+                              mfem::SparseMatrix* LORmat, double scale);
+   virtual ~AdditiveSchwarzLORSmoother() {}
+
+   virtual void SetOperator(const mfem::Operator& op) {}
+
+   virtual void Mult(const mfem::Vector& in, mfem::Vector& out) const;
+
+ private:
+   virtual void LocalSmoother(int element, const mfem::Vector& in,
+                              mfem::Vector& out) const;
+
+ protected:
+   const mfem::FiniteElementSpace& fespace_;
+   const Array<int>& ess_tdof_list_;
+   mfem::BilinearForm& aform_;
+   mfem::Vector diag_;
+   mfem::SparseMatrix* LORmat_;
+   double scale_;
+   mfem::Table el_to_el_;
+   Vector countingVector;
+   mutable DenseMatrixInverse inv;
+};
+
+
+class AdditiveSchwarzApproxLORSmoother : public Solver
+{
+ public:
+   AdditiveSchwarzApproxLORSmoother(const mfem::FiniteElementSpace& fespace,
+                              const Array<int>& ess_tdof_list,
+                              mfem::BilinearForm& aform,
+                              const mfem::Vector& diag,
+                              const mfem::Vector& LORdiag,
+                              mfem::SparseMatrix* LORmat, double scale);
+   virtual ~AdditiveSchwarzApproxLORSmoother() {}
+
+   virtual void SetOperator(const mfem::Operator& op) {}
+
+   virtual void Mult(const mfem::Vector& in, mfem::Vector& out) const;
+
+ private:
+   virtual void LocalSmoother(int element, const mfem::Vector& in,
+                              mfem::Vector& out) const;
+
+   void TensorProductMult2D(const DenseMatrix& A, const DenseMatrix& B, const Vector& in, Vector& out) const;   
+   void TensorProductMult2D(const Operator& A, const Operator& B, const Vector& in, Vector& out) const;
+
+ protected:
+   const mfem::FiniteElementSpace& fespace_;
+   const Array<int>& ess_tdof_list_;
+   mfem::BilinearForm& aform_;
+   mfem::Vector diag_;
+   mfem::Vector LORdiag_;
+   mfem::SparseMatrix* LORmat_;
+public:   
+   double scale_;
+   mfem::Table el_to_el_;
+   Vector countingVector;
+   DenseMatrix elmat;
+   Array<int> lexdofmap;
+   mutable std::vector<DenseMatrixInverse> inv;
+   std::vector<DenseMatrix> A1, A2, B1, B2;
+   std::vector<DenseMatrixInverse> invA2, invB1;
+   std::vector<DenseMatrixSchurDecomposition*> schurA1;
+   std::vector<DenseMatrixSchurDecomposition*> schurB2;
+   std::vector<DenseMatrixSylvesterSolver*> syl;
+};
+
+/**
+   Our first implementation is Jacobi within elements and Gauss-Seidel
+   between them.
+*/
+class ElementWiseJacobi : public ElementWiseSmoother
+{
+public:
+   ElementWiseJacobi(const mfem::FiniteElementSpace& fespace,
+                     mfem::BilinearForm& aform,
+                     const mfem::Vector& global_diag,
+                     double scale=1.0);
+
+   virtual void GetElementFromMatVec(int element, const mfem::Vector& x,
+                                     mfem::Vector& y) const;
+
+   virtual void LocalSmoother(int element, const mfem::Vector& in,
+                              mfem::Vector& out) const;
+
+protected:
+   mfem::BilinearForm& aform_;
+   const mfem::Vector& global_diag_;
+   double scale_;
+
+   mfem::Table el_to_el_;
+};
+
+}
+
+#endif // MFEM_ELEMENTWISESMOOTHER

fem/fem.hpp

@@ -33,6 +33,7 @@
 #include "tmop.hpp"
 #include "transfer.hpp"
 #include "spacehierarchy.hpp"
+#include "elementwisesmoother.hpp"
 #include "mgbilinearform.hpp"
 
 #ifdef MFEM_USE_MPI

fem/fespace.hpp

@@ -923,8 +923,8 @@ protected:
 
    mutable bool use_tensor_products;
 
-   static const int MAX_NQ2D = 100;
-   static const int MAX_ND2D = 100;
+   static const int MAX_NQ2D = 1000;
+   static const int MAX_ND2D = 1000;
    static const int MAX_VDIM2D = 2;
 
    static const int MAX_NQ3D = 1000;

general/forall.hpp

@@ -32,8 +32,8 @@ namespace mfem
 {
 
 // Maximum size of dofs and quads in 1D.
-const int MAX_D1D = 16;
-const int MAX_Q1D = 16;
+const int MAX_D1D = 32;
+const int MAX_Q1D = 32;
 
 // Implementation of MFEM's "parallel for" (forall) device/host kernel
 // interfaces supporting RAJA, CUDA, OpenMP, and sequential backends.

linalg/densemat.cpp

@@ -56,10 +56,23 @@ extern "C" void
 dgesvd_(char *JOBU, char *JOBVT, int *M, int *N, double *A, int *LDA,
         double *S, double *U, int *LDU, double *VT, int *LDVT, double *WORK,
         int *LWORK, int *INFO);
+extern "C" void
+dgehrd_(int *N, int *ILO, int *IHI, double *A, int *LDA, double* TAU,
+        double *WORK, int *LWORK, int *INFO);
+extern "C" void
+dorghr_(int *N, int *ILO, int *IHI, double *A, int *LDA, double* TAU,
+        double *WORK, int *LWORK, int *INFO);
+extern "C" void
+dhseqr_(char *JOB, char *COMPZ, int *N, int *ILO, int *IHI, double *H,
+        int *LDH, double *WR, double *WI, double *Z, int *LDZ, double *WORK,
+        int *LWORK, int *INFO);
+extern "C" void 
+dtrsyl_(char* TRANA, char* TRANB, int* ISGN, int* M, int* N,
+        double* A, int* LDA, double* B, int* LDB, double* C,
+        int* LDC, double* SCALE, int* INFO);
 #endif
 
-
-namespace mfem
+    namespace mfem
 {
 
 using namespace std;
@@ -4403,7 +4416,9 @@ void DenseMatrixSVD::Init()
 {
 #ifdef MFEM_USE_LAPACK
    sv.SetSize(min(m, n));
-
+   U.SetSize(m,m);
+   V_T.SetSize(n,n);
+   
    jobu  = 'N';
    jobvt = 'N';
 
@@ -4429,8 +4444,10 @@ void DenseMatrixSVD::Eval(DenseMatrix &M)
 #endif
 
 #ifdef MFEM_USE_LAPACK
-   dgesvd_(&jobu, &jobvt, &m, &n, M.Data(), &m, sv.GetData(), NULL, &m,
-           NULL, &n, work, &lwork, &info);
+   jobu  = 'A';
+   jobvt = 'A';
+   dgesvd_(&jobu, &jobvt, &m, &n, M.Data(), &m, sv.GetData(), U.Data(), &m,
+           V_T.Data(), &n, work, &lwork, &info);
 
    if (info)
    {
@@ -4449,6 +4466,108 @@ DenseMatrixSVD::~DenseMatrixSVD()
 #endif
 }
 
+DenseMatrixSchurDecomposition::DenseMatrixSchurDecomposition(const DenseMatrix &M)
+{
+#ifdef MFEM_USE_LAPACK
+   int n = M.Height();
+   T = M;
+   int ILO = 1;
+   int IHI = n;
+   int info = 0;
+
+   Vector tau(n - 1);
+
+   double work_;
+   int lwork = -1;
+   dgehrd_(&n, &ILO, &IHI, T.Data(), &n, tau.GetData(), &work_, &lwork, &info);
+
+   lwork = static_cast<int>(work_);
+   double* work = new double[lwork];
+   dgehrd_(&n, &ILO, &IHI, T.Data(), &n, tau.GetData(), work, &lwork, &info);
+   delete[] work;
+
+   if (info)
+   {
+      mfem::err << "dgehrd_ : info = " << info << endl;
+      mfem_error();
+   }
+
+   Q = T;
+   lwork = -1;
+   dorghr_(&n, &ILO, &IHI, Q.Data(), &n, tau.GetData(), &work_, &lwork,
+           &info);
+   lwork = static_cast<int>(work_);
+   work = new double[lwork];
+   dorghr_(&n, &ILO, &IHI, Q.Data(), &n, tau.GetData(), work, &lwork, &info);
+   delete[] work;
+
+   if (info)
+   {
+      mfem::err << "dorghr_ : info = " << info << endl;
+      mfem_error();
+   }
+
+   char job = 'S';
+   char compz = 'V';
+   Vector wr(n);
+   Vector wi(n);
+   lwork = -1;
+   dhseqr_(&job, &compz, &n, &ILO, &IHI, T.Data(), &n, wr.GetData(),
+           wi.GetData(), Q.Data(), &n, &work_, &lwork, &info);
+   lwork = static_cast<int>(work_);
+   work = new double[lwork];
+   dhseqr_(&job, &compz, &n, &ILO, &IHI, T.Data(), &n, wr.GetData(),
+           wi.GetData(), Q.Data(), &n, work, &lwork, &info);
+   delete[] work;
+
+   if (info)
+   {
+      mfem::err << "dhseqr_ : info = " << info << endl;
+      mfem_error();
+   }
+
+   Q_T = Q;
+   Q_T.Transpose();
+
+#else
+   mfem_error("DenseMatrixSchurDecomposition::Eval(): Compiled without LAPACK");
+#endif
+}
+
+DenseMatrixSylvesterSolver::DenseMatrixSylvesterSolver(const DenseMatrix& A_, const DenseMatrix& B_, bool tranA_, bool tranB_)
+   : A(A_), B(B_)
+{
+   tranA = tranA_ ? 'T' : 'N';
+   tranB = tranB_ ? 'T' : 'N';
+}
+
+void DenseMatrixSylvesterSolver::Mult(const Vector &c, Vector &x) const
+{
+   int isgn = 1;
+   int m = A.Width();
+   int n = B.Height();
+   double scale = 1.0;
+   int info = 0;
+
+   x = c;
+   dtrsyl_(&tranA, &tranB, &isgn, &m, &n, A.Data(), &m, B.Data(), &n, x.GetData(), &m, &scale, &info);
+
+   if (info)
+   {
+      mfem::err << "dtrsyl_ : info = " << info << endl;
+      mfem_error();
+   }
+
+   // Test inversion
+   // DenseMatrix C(c.GetData(), m, n);
+   // DenseMatrix X(x.GetData(), m, n);
+   // DenseMatrix tmp(m,n);
+   // mfem::Mult(A, X, tmp);
+   // mfem::AddMultABt(X, B, tmp);
+   // tmp -= C;
+   // std::cout << "test = " << tmp.MaxMaxNorm() << std::endl;
+}
+
 
 void DenseTensor::AddMult(const Table &elem_dof, const Vector &x, Vector &y)
 const

linalg/densemat.hpp

@@ -631,6 +631,8 @@ public:
 class DenseMatrixSVD
 {
    Vector sv;
+   DenseMatrix U;
+   DenseMatrix V_T;
    int m, n;
 
 #ifdef MFEM_USE_LAPACK
@@ -647,9 +649,48 @@ public:
    void Eval(DenseMatrix &M);
    Vector &Singularvalues() { return sv; }
    double Singularvalue(int i) { return sv(i); }
+   DenseMatrix& GetU() { return U; }
+   DenseMatrix& GetV_T() { return V_T; }
+
    ~DenseMatrixSVD();
 };
 
+class DenseMatrixSchurDecomposition
+{
+   DenseMatrix T;
+   DenseMatrix Q;
+   DenseMatrix Q_T;
+
+public:
+   DenseMatrixSchurDecomposition(const DenseMatrix &M);
+   DenseMatrix& GetT() { return T; };
+   DenseMatrix& GetQ() { return Q; };
+   DenseMatrix& GetQ_T() { return Q_T; };
+};
+
+/// Solver for the matrix equation op(A)*X + X*op(B) = C where
+/// op(A) = A or A**T using Bartels–Stewart algorithm. A and B must be
+/// quasi-triangular matrices.
+class DenseMatrixSylvesterSolver
+{
+ private:
+   const DenseMatrix& A;
+   const DenseMatrix& B;
+
+#ifdef MFEM_USE_LAPACK
+   mutable char tranA;
+   mutable char tranB;
+#endif
+
+ public:
+   DenseMatrixSylvesterSolver(const DenseMatrix& A, const DenseMatrix& B,
+                              bool tranA_ = false, bool tranB_ = false);
+
+   /// c = vec(C)
+   /// x = vec(X)
+   void Mult(const Vector& c, Vector& x) const;
+};
+
 class Table;
 
 /// Rank 3 tensor (array of matrices)

linalg/solvers.cpp

@@ -124,7 +124,7 @@ void OperatorJacobiSmoother::Setup()
    auto X = dinv.Write();
    MFEM_FORALL(i, N, X[i] = delta / D[i]; );
    auto I = ess_tdof_list.Read();
-   MFEM_FORALL(i, ess_tdof_list.Size(), X[I[i]] = delta; );
+   MFEM_FORALL(i, ess_tdof_list.Size(), X[I[i]] = 1.0; );
 }
 
 void OperatorJacobiSmoother::Mult(const Vector& x, Vector &y) const

miniapps/performance/mg.cpp

@@ -1,10 +1,18 @@
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
+#include "../../general/forall.hpp"
 
 using namespace std;
 using namespace mfem;
 
+static const double omega = 8.0 * M_PI;
+
+double exactFun(Vector& x)
+{
+   return std::sin(omega * (x[0]+x[1]) / std::sqrt(2));
+}
+
 class PoissonMultigridOperator : public TimedMultigridOperator
 {
  private:
@@ -19,11 +27,13 @@ class PoissonMultigridOperator : public TimedMultigridOperator
    ParFiniteElementSpace* fespace_lor{nullptr};
    bool useCoarsePCG{false};
    HypreBoomerAMG* amg{nullptr};
-   CGSolver* coarsePCGSolver{nullptr};
+   GMRESSolver* coarsePCGSolver{nullptr};
 
-   void AddIntegrators(ParBilinearForm* form)
+   void AddIntegrators(BilinearForm* form)
    {
       form->AddDomainIntegrator(new DiffusionIntegrator(*coeff));
+      ConstantCoefficient* massCoeff = new ConstantCoefficient(-omega * omega);
+      form->AddDomainIntegrator(new MassIntegrator(*massCoeff));
    }
 
    Operator* ConstructOperator(ParFiniteElementSpace* fespace,
@@ -100,22 +110,22 @@ class PoissonMultigridOperator : public TimedMultigridOperator
       amg->SetPrintLevel(-1);
       amg->SetMaxIter(coarseSteps);
 
-      if (useCoarsePCG)
-      {
+      // if (useCoarsePCG)
+      // {
          amg->SetMaxIter(1);
-         coarsePCGSolver = new CGSolver(MPI_COMM_WORLD);
-         coarsePCGSolver->SetPrintLevel(-1);
-         coarsePCGSolver->SetMaxIter(50);
-         coarsePCGSolver->SetRelTol(1e-3);
+         coarsePCGSolver = new GMRESSolver(MPI_COMM_WORLD);
+         coarsePCGSolver->SetPrintLevel(0);
+         coarsePCGSolver->SetMaxIter(5000);
+         coarsePCGSolver->SetRelTol(1e-12);
          coarsePCGSolver->SetAbsTol(0.0);
          coarsePCGSolver->SetOperator(*opr);
-         coarsePCGSolver->SetPreconditioner(*amg);
+         // coarsePCGSolver->SetPreconditioner(*amg);
          return coarsePCGSolver;
-      }
-      else
-      {
-         return amg;
-      }
+      // }
+      // else
+      // {
+      //    return amg;
+      // }
    }
 
  public:
@@ -129,7 +139,8 @@ class PoissonMultigridOperator : public TimedMultigridOperator
       if (jump)
       {
          auto f = [](const Vector& x) {
-            return (x[0] + x[1] < 1.0) ? 1.0 : 100.0;
+            // return (x[0] < 0.5) ? 1.0 : 100.0;
+            return 5 * x[0] + 1.0;
          };
          coeff = new FunctionCoefficient(f);
       }
@@ -178,20 +189,96 @@ class PoissonMultigridOperator : public TimedMultigridOperator
 
       if (partialAssembly)
       {
-         Vector diag(fespace->GetTrueVSize());
-         forms.Last()->AssembleDiagonal(diag);
+         Vector* diag = new Vector(fespace->GetTrueVSize());
+         forms.Last()->AssembleDiagonal(*diag);
 
-         Vector ev(solveOperator->Width());
-         OperatorJacobiSmoother invDiagOperator(diag, essentialDofs, 1.0);
-         ProductOperator diagPrecond(&invDiagOperator, solveOperator, false,
-                                     false);
+         Vector* coeffDiag = new Vector(fespace->GetTrueVSize());
+         {
+            Array<int> local_dofs;
+            int ne = fespace->GetNE();
+            const IntegrationRule& ir = fespace->GetFE(0)->GetNodes();
+            int nq = ir.GetNPoints();
+
+            for(int e = 0; e < ne; ++e)
+            {
+               fespace->GetElementDofs(e, local_dofs);
+               ElementTransformation& T = *fespace->GetElementTransformation(e);
+               for (int q = 0; q < nq; ++q)
+               {
+                  (*coeffDiag)[local_dofs[q]] = 1.0 / sqrt(coeff->Eval(T, ir.IntPoint(q)));
+               }
+            }
+         }
+
+         Mesh* pmesh_lor = new Mesh(fespace->GetMesh(), fespace->GetOrder(0), BasisType::GaussLobatto);
+         H1_FECollection* fec_lor_ = new H1_FECollection(1, pmesh_lor->Dimension(), BasisType::GaussLobatto);
+         FiniteElementSpace* fespace_lor_ = new FiniteElementSpace(pmesh_lor, fec_lor_);
+
+         BilinearForm* a_pc_ = new BilinearForm(fespace_lor_);
+         a_pc_->SetAssemblyLevel(AssemblyLevel::FULL);
+         AddIntegrators(a_pc_);
+         // ConstantCoefficient* constCoeff = new ConstantCoefficient(1.0);
+         // a_pc_->AddDomainIntegrator(new DiffusionIntegrator(*constCoeff));
+         a_pc_->UsePrecomputedSparsity();
+         a_pc_->Assemble();
+         
+         // BilinearForm* a_pc_pa = new BilinearForm(fespace_lor_);
+         // a_pc_pa->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+         // AddIntegrators(a_pc_pa);
+         // a_pc_pa->Assemble();
+         Vector* LORdiag = new Vector(fespace->GetTrueVSize());
+         // a_pc_pa->AssembleDiagonal(*LORdiag);
+
+         SparseMatrix* LORmat = new SparseMatrix();
+         a_pc_->FormSystemMatrix(essentialDofs, *LORmat);
 
          PowerMethod powerMethod(MPI_COMM_WORLD);
-         double estLargestEigenvalue =
-             powerMethod.EstimateLargestEigenvalue(diagPrecond, ev, 10, 1e-8);
-         smoother = new OperatorChebyshevSmoother(solveOperator, diag,
-                                                  essentialDofs, chebyshevOrder,
-                                                  estLargestEigenvalue);
+
+         // AdditiveSchwarzApproxLORSmoother test(*fespace, essentialDofs, *forms.Last(), *coeffDiag, LORmat, 1.0);
+         AdditiveSchwarzApproxLORSmoother& test = *new AdditiveSchwarzApproxLORSmoother(*fespace, essentialDofs, *forms.Last(), *coeffDiag, *LORdiag, LORmat, 1.0);
+         // OperatorJacobiSmoother test(*diag, essentialDofs, 1.0);
+
+
+         ProductOperator powerOperator(&test, solveOperator, false, false);
+         Vector ev(solveOperator->Width());
+         double estLargestEigenvalue = powerMethod.EstimateLargestEigenvalue(powerOperator, ev, 10, 1e-8);
+
+         std::cout << "ev = " << estLargestEigenvalue << std::endl;
+
+         double upper_bound = 1.1 * estLargestEigenvalue;
+         double lower_bound = 0.0 * estLargestEigenvalue;
+         double theta = 0.5 * (upper_bound + lower_bound);
+         double delta = 0.5 * (upper_bound - lower_bound);
+         double weight = 1.0 / theta;
+         std::cout << "weight = " << weight << std::endl;
+
+         // HypreParMatrix* LORmatp = new HypreParMatrix();
+         // a_pc_->FormSystemMatrix(essentialDofs, *LORmatp);
+
+         // SparseMatrix* LORmat = new SparseMatrix();
+         // LORmatp->GetDiag(*LORmat);
+
+         std::cout << "truevsize = " << fespace->GetTrueVSize() << std::endl;
+         std::cout << "Width = " << LORmat->Width() << std::endl;
+
+         smoother = &test; test.scale_ = weight;
+         // smoother = new AdditiveSchwarzApproxLORSmoother(*fespace, essentialDofs, *forms.Last(), *coeffDiag, *LORdiag, LORmat, weight);
+         // smoother = new ElementWiseJacobi(*fespace, *forms.Last(), *diag, essentialDofs, 2.0/3.0);
+         // smoother = new OperatorJacobiSmoother(*diag, essentialDofs, weight);
+         // smoother = new OperatorChebyshevSmoother(solveOperator, *diag, essentialDofs, chebyshevOrder, estLargestEigenvalue);
+         
+
+         // Vector ev(solveOperator->Width());
+         // OperatorJacobiSmoother invDiagOperator(*diag, essentialDofs, 1.0);
+         // ProductOperator diagPrecond(&invDiagOperator, solveOperator, false,
+         //                             false);
+
+         // PowerMethod powerMethod(MPI_COMM_WORLD);
+         // double estLargestEigenvalue =
+         //     powerMethod.EstimateLargestEigenvalue(diagPrecond, ev, 10, 1e-8);
+         // smoother = new OperatorChebyshevSmoother(solveOperator, *diag,
+         //                                          essentialDofs, chebyshevOrder,
+         //                                          estLargestEigenvalue);
       }
       else
       {
@@ -236,11 +323,11 @@ int main(int argc, char* argv[])
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
 
    // 1. Parse command-line options.
-   const char* mesh_file = "../../data/inline-hex.mesh";
+   const char* mesh_file = "../../data/inline-quad.mesh";
    int ref_levels = 0;
    int pref_levels = 0;
    int order = 1;
-   int h_levels = 1;
+   int h_levels = 2;
    int o_levels = 1;
    int smoothingSteps = 3;
    int coarseSteps = 2;
@@ -339,7 +426,8 @@ int main(int argc, char* argv[])
    // 2. Read the mesh from the given mesh file. We can handle triangular,
    //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
    //    the same code.
-   Mesh* mesh = new Mesh(mesh_file, 1, 1);
+   // Mesh* mesh = new Mesh(mesh_file, 1, 1);
+   Mesh* mesh = new Mesh(1,1, Element::QUADRILATERAL, true, 1.0, 1.0, false);
    int dim = mesh->Dimension();
 
    Array<int> ess_bdr(mesh->bdr_attributes.Max());
@@ -361,6 +449,8 @@ int main(int argc, char* argv[])
       pmesh->UniformRefinement();
    }
 
+   std::cout << "Number of elements: " << pmesh->GetNE() << std::endl;
+
    Array<int> orders;
    Array<FiniteElementCollection*> feCollectons;
    orders.Append(order);
@@ -495,6 +585,9 @@ int main(int argc, char* argv[])
    ParGridFunction x(&spaceHierarchy->GetFinestFESpace());
    x = 0.0;
 
+   FunctionCoefficient exact(exactFun);
+   x.ProjectCoefficient(exact);
+
    if (myid == 0)
    {
       cout << "Assembling rhs..." << flush;
@@ -503,7 +596,7 @@ int main(int argc, char* argv[])
    tic_toc.Start();
    ParLinearForm* b = new ParLinearForm(&spaceHierarchy->GetFinestFESpace());
    ConstantCoefficient one(1.0);
-   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   // b->AddDomainIntegrator(new DomainLFIntegrator(one));
    b->Assemble();
    tic_toc.Stop();
    if (myid == 0)
@@ -517,13 +610,15 @@ int main(int argc, char* argv[])
    tic_toc.Clear();
    tic_toc.Start();
 
-   CGSolver pcg(MPI_COMM_WORLD);
+   GMRESSolver pcg(MPI_COMM_WORLD);
+   pcg.SetKDim(100);
    pcg.SetPrintLevel(1);
-   pcg.SetMaxIter(100);
-   pcg.SetRelTol(1e-6);
+   pcg.SetMaxIter(1000);
+   pcg.SetRelTol(1e-10);
    pcg.SetAbsTol(0.0);
    pcg.SetOperator(*solveOperator);
    pcg.SetPreconditioner(*preconditioner);
+   // pcg.SetPreconditioner(*solveOperator->GetSmootherAtLevel(spaceHierarchy->GetFinestLevelIndex()));
    pcg.Mult(B, X);
 
    tic_toc.Stop();
@@ -536,18 +631,15 @@ int main(int argc, char* argv[])
       if (TimedMultigridOperator* tmg =
               dynamic_cast<TimedMultigridOperator*>(solveOperator))
       {
-         tmg->PrintStats(TimedMultigridOperator::Operation::OPERATOR, cout);
-         tmg->PrintStats(TimedMultigridOperator::Operation::PROLONGATION, cout);
-         tmg->PrintStats(TimedMultigridOperator::Operation::RESTRICTION, cout);
-         tmg->PrintStats(TimedMultigridOperator::Operation::SMOOTHER, cout);
+         // tmg->PrintStats(TimedMultigridOperator::Operation::OPERATOR, cout);
+         // tmg->PrintStats(TimedMultigridOperator::Operation::PROLONGATION, cout);
+         // tmg->PrintStats(TimedMultigridOperator::Operation::RESTRICTION, cout);
+         // tmg->PrintStats(TimedMultigridOperator::Operation::SMOOTHER, cout);
       }
    }
 
    solveOperator->RecoverFEMSolution(X, *b, x);
 
-   delete preconditioner;
-   delete solveOperator;
-
    if (visualization)
    {
       char vishost[] = "localhost";
@@ -560,6 +652,8 @@ int main(int argc, char* argv[])
                << flush;
    }
 
+   delete preconditioner;
+   delete solveOperator;
    delete b;
    delete spaceHierarchy;
 

tests/unit/fem/test_elementmatmult.cpp

@@ -0,0 +1,105 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+using namespace mfem;
+
+namespace elementmatmult
+{
+
+TEST_CASE("elementmatmult")
+{
+   for (int dimension = 2; dimension < 4; ++dimension)
+   {
+      for (int integrator = 0; integrator < 3; ++integrator)
+      {
+         for (int ne = 1; ne < 3; ++ne)
+         {
+            std::cout << "Testing " << dimension << "D partial assembly element mat mult: "
+                      << std::pow(ne, dimension) << " elements,"
+                      << " integrator " << integrator << std::endl;
+            for (int order = 1; order < 5; ++order)
+            {
+               Mesh * mesh;
+               if (dimension == 2)
+               {
+                  mesh = new Mesh(ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
+               }
+               else
+               {
+                  mesh = new Mesh(ne, ne, ne, Element::HEXAHEDRON, 1, 1.0, 1.0, 1.0);
+               }
+               FiniteElementCollection *h1_fec = new H1_FECollection(order, dimension);
+               FiniteElementSpace h1_fespace(mesh, h1_fec);
+               Array<int> ess_tdof_list;
+               // Array<int> ess_bdr(mesh->bdr_attributes.Max());
+               // ess_bdr = 1;
+               // h1_fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+               BilinearForm paform(&h1_fespace);
+               ConstantCoefficient one(1.0);
+               paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+               if (integrator < 2)
+               {
+                  paform.AddDomainIntegrator(new DiffusionIntegrator(one));
+               }
+               if (integrator > 0)
+               {
+                  paform.AddDomainIntegrator(new MassIntegrator(one));
+               }
+               paform.Assemble();
+
+               GridFunction x(&h1_fespace);
+               x = 0.0;
+               GridFunction b(&h1_fespace);
+               b = 1.0;
+               BilinearForm assemblyform(&h1_fespace);
+               if (integrator < 2)
+               {
+                  assemblyform.AddDomainIntegrator(new DiffusionIntegrator(one));
+               }
+               if (integrator > 0)
+               {
+                  assemblyform.AddDomainIntegrator(new MassIntegrator(one));
+               }
+               // assemblyform.SetDiagonalPolicy(Matrix::DIAG_ONE); // todo?
+               assemblyform.Assemble();
+               assemblyform.Finalize();
+               OperatorPtr A_assembly;
+               Vector B, X;
+               assemblyform.FormLinearSystem(ess_tdof_list, x, b, A_assembly, X, B);
+
+               for (int elno = 0; elno < std::min(8, mesh->GetNE()); ++elno)
+               {
+                  Vector xin(h1_fespace.GetTrueVSize());
+                  xin.Randomize();
+                  Vector y_assembly, y_pa;
+                  assemblyform.ElementMatrixMult(elno, xin, y_assembly);
+                  paform.ElementMatrixMult(elno, xin, y_pa);
+
+                  y_assembly -= y_pa;
+                  double error = y_assembly.Norml2();
+                  std::cout << "    order: " << order << ", elno " << elno
+                            << ", error norm: " << error << std::endl;
+                  REQUIRE(y_assembly.Norml2() < 1.e-14);
+               }
+
+               delete mesh;
+               delete h1_fec;
+            }
+         }
+      }
+   }
+}
+
+} // namespace elementmatmult

tests/unit/fem/test_elementwisesmoother.cpp

@@ -0,0 +1,109 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+using namespace mfem;
+
+namespace elementwisesmoother
+{
+
+TEST_CASE("elementwisesmoother")
+{
+   for (int dimension = 2; dimension < 4; ++dimension)
+   {
+      for (int integrator = 0; integrator < 2; ++integrator)
+      {
+         const int ne = 1;
+         std::cout << "Testing " << dimension << "D element-wise smoother with "
+                   << "integrator " << integrator << std::endl;
+         for (int order = 1; order < 5; ++order)
+         {
+            Mesh * mesh;
+            if (dimension == 2)
+            {
+               mesh = new Mesh(ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
+            }
+            else
+            {
+               mesh = new Mesh(ne, ne, ne, Element::HEXAHEDRON, 1, 1.0, 1.0, 1.0);
+            }
+            FiniteElementCollection *h1_fec = new H1_FECollection(order, dimension);
+            FiniteElementSpace h1_fespace(mesh, h1_fec);
+            Array<int> ess_tdof_list;
+ 
+            BilinearForm paform(&h1_fespace);
+            ConstantCoefficient one(1.0); 
+            paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+            if (integrator < 2)
+            {
+               paform.AddDomainIntegrator(new DiffusionIntegrator(one));
+            }
+            if (integrator > 0)
+            {
+               paform.AddDomainIntegrator(new MassIntegrator(one));
+            }
+            paform.Assemble();
+
+            BilinearForm assemblyform(&h1_fespace);
+            if (integrator < 2)
+            {
+               assemblyform.AddDomainIntegrator(new DiffusionIntegrator(one));
+            }
+            if (integrator > 0)
+            {
+               assemblyform.AddDomainIntegrator(new MassIntegrator(one));
+            }
+            assemblyform.SetDiagonalPolicy(Matrix::DIAG_ONE);
+            assemblyform.Assemble();
+            assemblyform.Finalize();
+            const SparseMatrix& A_explicit = assemblyform.SpMat();
+
+            Vector pa_diag(h1_fespace.GetVSize());
+            paform.AssembleDiagonal(pa_diag);
+            ElementWiseJacobi e_pa_jacobi(h1_fespace, paform, pa_diag);
+            ElementWiseJacobi e_assembly_jacobi(h1_fespace, assemblyform, pa_diag);
+            DSmoother mat_jacobi(A_explicit);
+
+            Vector xin(h1_fespace.GetTrueVSize());
+            xin.Randomize();
+            Vector y_mat(xin);
+            y_mat = 0.0;
+            Vector y_assembly(xin);
+            y_assembly = 0.0;
+            Vector y_pa(xin);
+            y_pa = 0.0;
+
+            e_pa_jacobi.Mult(xin, y_pa);
+            e_assembly_jacobi.Mult(xin, y_assembly);
+            mat_jacobi.Mult(xin, y_mat);
+
+            y_pa -= y_mat;
+            double pa_error = y_pa.Norml2();
+            std::cout << "  order: " << order << ", pa error norm: " << pa_error
+                      << std::endl;
+            REQUIRE(pa_error < 1.e-12);
+
+            y_assembly -= y_mat;
+            double assembly_error = y_assembly.Norml2();
+            std::cout << "  order: " << order << ", assembly error norm: "
+                      << assembly_error << std::endl;
+            REQUIRE(assembly_error < 1.e-12);
+
+            delete mesh;
+            delete h1_fec;
+         }
+      }
+   }
+}
+
+} // namespace elementwisesmoother