From 6b5db0c503c8899930329332769f12011006ed4c Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Tue, 12 Sep 2023 11:55:15 -0700 Subject: [PATCH] Add linalg kernels::CalcLeftInverse --- fem/qinterp/grad.hpp | 33 ++++++++++----------------------- linalg/densemat.cpp | 27 ++++----------------------- linalg/kernels.hpp | 39 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 53 insertions(+), 46 deletions(-) diff --git a/fem/qinterp/grad.hpp b/fem/qinterp/grad.hpp index 61228b8aace..8f0a2d3d2a0 100644 --- a/fem/qinterp/grad.hpp +++ b/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; diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp index 15ef7c0cfc1..86391f2397f 100644 --- a/linalg/densemat.cpp +++ b/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 diff --git a/linalg/kernels.hpp b/linalg/kernels.hpp index 1c05e84b16f..62230548e02 100644 --- a/linalg/kernels.hpp +++ b/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 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.