Add linalg kernels::CalcLeftInverse

fem/qinterp/grad.hpp

@@ -62,19 +62,21 @@ static void Derivatives1D(const int NE,
  else if (sdim == 2)
  {
  const double Jloc[2] = {j(q,0,e), j(q,1,e)};
- const double t = 1.0 / (Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1]);
- const double U = t*Jloc[0]*du[0];
- const double V = t*Jloc[1]*du[0];
+ double Jinv[3];
+ kernels::CalcLeftInverse<2,1>(Jloc, Jinv);
+ const double U = Jinv[0]*du[0];
+ const double V = Jinv[1]*du[0];
  du[0] = U;
  du[1] = V;
  }
  else // sdim == 3
  {
  const double Jloc[3] = {j(q,0,e), j(q,1,e), j(q,2,e)};
- const double t = 1.0/(Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1] + Jloc[2]*Jloc[2]);
- const double U = t*Jloc[0]*du[0];
- const double V = t*Jloc[1]*du[0];
- const double W = t*Jloc[2]*du[0];
+ double Jinv[3];
+ kernels::CalcLeftInverse<3,1>(Jloc, Jinv);
+ const double U = Jinv[0]*du[0];
+ const double V = Jinv[1]*du[0];
+ const double W = Jinv[2]*du[0];
  du[0] = U;
  du[1] = V;
  du[2] = W;
@@ -194,25 +196,10 @@ static void Derivatives2D(const int NE,
  Jloc[3] = j(qx,qy,0,1,e);
  Jloc[4] = j(qx,qy,1,1,e);
  Jloc[5] = j(qx,qy,2,1,e);
-
- double ee, gg, ff;
- ee = Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1] + Jloc[2]*Jloc[2];
- gg = Jloc[3]*Jloc[3] + Jloc[4]*Jloc[4] + Jloc[5]*Jloc[5];
- ff = Jloc[0]*Jloc[3] + Jloc[1]*Jloc[4] + Jloc[2]*Jloc[5];
- const double t = 1.0 / (ee*gg - ff*ff);
- ee *= t; gg *= t; ff *= t;
-
- Jinv[0] = Jloc[0]*gg - Jloc[3]*ff;
- Jinv[1] = Jloc[3]*ee - Jloc[0]*ff;
- Jinv[2] = Jloc[1]*gg - Jloc[4]*ff;
- Jinv[3] = Jloc[4]*ee - Jloc[1]*ff;
- Jinv[4] = Jloc[2]*gg - Jloc[5]*ff;
- Jinv[5] = Jloc[5]*ee - Jloc[2]*ff;
-
+ kernels::CalcLeftInverse<3,2>(Jloc, Jinv);
  const double U = Jinv[0]*du[0] + Jinv[1]*du[1];
  const double V = Jinv[2]*du[0] + Jinv[3]*du[1];
  const double W = Jinv[4]*du[0] + Jinv[5]*du[1];
-
  du[0] = U;
  du[1] = V;
  du[2] = W;

linalg/densemat.cpp

@@ -2589,49 +2589,30 @@ void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
  MFEM_ASSERT(inva.Height() == a.Width(), "incorrect dimensions");
  MFEM_ASSERT(inva.Width() == a.Height(), "incorrect dimensions");
 
- double t;
-
  if (a.Width() < a.Height())
  {
  const double *d = a.Data();
  double *id = inva.Data();
  if (a.Height() == 2)
  {
- t = 1.0 / (d[0]*d[0] + d[1]*d[1]);
- id[0] = d[0] * t;
- id[1] = d[1] * t;
+ kernels::CalcLeftInverse<2,1>(d, id);
  }
  else
  {
  if (a.Width() == 1)
  {
- t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
- id[0] = d[0] * t;
- id[1] = d[1] * t;
- id[2] = d[2] * t;
+ kernels::CalcLeftInverse<3,1>(d, id);
  }
  else
  {
- double e, g, f;
- e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
- g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
- f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
- t = 1.0 / (e*g - f*f);
- e *= t; g *= t; f *= t;
-
- id[0] = d[0]*g - d[3]*f;
- id[1] = d[3]*e - d[0]*f;
- id[2] = d[1]*g - d[4]*f;
- id[3] = d[4]*e - d[1]*f;
- id[4] = d[2]*g - d[5]*f;
- id[5] = d[5]*e - d[2]*f;
+ kernels::CalcLeftInverse<3,2>(d, id);
  }
  }
  return;
  }
 
 #ifdef MFEM_DEBUG
- t = a.Det();
+ const double t = a.Det();
  MFEM_ASSERT(std::abs(t) > 1.0e-14 * pow(a.FNorm()/a.Width(), a.Width()),
  "singular matrix!");
 #endif

linalg/kernels.hpp

@@ -406,6 +406,10 @@ void MultAtB(const int Aheight, const int Awidth, const int Bwidth,
  }
 }
 
+/// Given a matrix of size 2x1, 3x1, or 3x2, compute the left inverse.
+template<int HEIGHT, int WIDTH> MFEM_HOST_DEVICE
+void CalcLeftInverse(const double *data, double *left_inv);
+
 /// Compute the spectrum of the matrix of size dim with given @a data, returning
 /// the eigenvalues in the array @a lambda and the eigenvectors in the array @a
 /// vec (listed consecutively).
@@ -1062,6 +1066,41 @@ int Reduce3S(const int &mode,
 
 } // namespace kernels::internal
 
+// Implementations of CalcLeftInverse for dim = 1, 2.
+
+template<> MFEM_HOST_DEVICE inline
+void CalcLeftInverse<2,1>(const double *d, double *left_inv)
+{
+ const double t = 1.0 / (d[0]*d[0] + d[1]*d[1]);
+ left_inv[0] = d[0] * t;
+ left_inv[1] = d[1] * t;
+}
+
+template<> MFEM_HOST_DEVICE inline
+void CalcLeftInverse<3,1>(const double *d, double *left_inv)
+{
+ const double t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
+ left_inv[0] = d[0] * t;
+ left_inv[1] = d[1] * t;
+ left_inv[2] = d[2] * t;
+}
+
+template<> MFEM_HOST_DEVICE inline
+void CalcLeftInverse<3,2>(const double *d, double *left_inv)
+{
+ double e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+ double g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+ double f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+ const double t = 1.0 / (e*g - f*f);
+ e *= t; g *= t; f *= t;
+
+ left_inv[0] = d[0]*g - d[3]*f;
+ left_inv[1] = d[3]*e - d[0]*f;
+ left_inv[2] = d[1]*g - d[4]*f;
+ left_inv[3] = d[4]*e - d[1]*f;
+ left_inv[4] = d[2]*g - d[5]*f;
+ left_inv[5] = d[5]*e - d[2]*f;
+}
 
 // Implementations of CalcEigenvalues and CalcSingularvalue for dim = 2, 3.