[!295] Feature/loopy invert matrix

Merge branch 'feature/loopy-invert-matrix' into 'master'

ref:extensions/dune-codegen Moves the matrix inversion stuff from
blockstructured/geometry.py to pdelab/tensor.py. This touches on the subject
of [#139].

You can get the name of the inverse of a matrix A with shape shape from

    name = name_matrix_inverse(A, shape, visitor)

matrix_inverse() returns the inverse wrapped in a pymbolic variable, although
I don't know how useful this is.

TODO:

-   [x] handle inames correctly
-   [ ] what about tags? -\> not really necessary
-   [ ] what about dependencies -\> not really necessary

See merge request [extensions/dune-codegen!295]

  [#139]: gitlab.dune-project.org/NoneNone/issues/139
  [extensions/dune-codegen!295]: gitlab.dune-project.org/extensions/dune-codegen/merge_requests/295

python/dune/codegen/blockstructured/basis.py

@@ -22,7 +22,7 @@ from dune.codegen.pdelab.geometry import world_dimension, component_iname
 from dune.codegen.pdelab.spaces import type_leaf_gfs, name_lfs
 from dune.codegen.pdelab.restriction import restricted_name
 from dune.codegen.blockstructured.spaces import lfs_inames
-from dune.codegen.blockstructured.tools import tensor_index_to_sequential_index, sub_element_inames
+from dune.codegen.blockstructured.tools import tensor_index_to_sequential_index, remove_sub_element_inames
 
 from ufl import MixedElement
 
@@ -64,7 +64,7 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
         cache = name_localbasis_cache(element)
         qp = self.to_cell(self.quadrature_position_in_micro())
         localbasis = name_localbasis(element)
-        instruction(inames=self.quadrature_inames(),
+        instruction(inames=remove_sub_element_inames(self.quadrature_inames()),
                     code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
                     assignees=frozenset({name}),
                     read_variables=frozenset({get_pymbolic_basename(qp)}),
@@ -100,7 +100,7 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
         cache = name_localbasis_cache(element)
         qp = self.to_cell(self.quadrature_position_in_micro())
         localbasis = name_localbasis(element)
-        instruction(inames=self.quadrature_inames(),
+        instruction(inames=remove_sub_element_inames(self.quadrature_inames()),
                     code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
                     assignees=frozenset({name}),
                     read_variables=frozenset({get_pymbolic_basename(qp)}),
@@ -132,7 +132,7 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
 
         instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
                     assignee=assignee,
-                    forced_iname_deps=frozenset(self.quadrature_inames() + (dimindex,) + sub_element_inames()),
+                    forced_iname_deps=frozenset(self.quadrature_inames() + (dimindex,)),
                     forced_iname_deps_is_final=True,
                     )
 
@@ -159,7 +159,7 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
 
         instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
                     assignee=assignee,
-                    forced_iname_deps=frozenset(self.quadrature_inames() + sub_element_inames()),
+                    forced_iname_deps=frozenset(self.quadrature_inames()),
                     forced_iname_deps_is_final=True,
                     )
 

python/dune/codegen/blockstructured/geometry.py

+import pymbolic.primitives as prim
+from loopy.match import Writes
+
+from dune.codegen.blockstructured.tools import name_point_in_macro
 from dune.codegen.generation import (geometry_mixin,
                                      temporary_variable,
                                      instruction,
                                      get_global_context_value,
                                      domain)
-from dune.codegen.tools import get_pymbolic_basename
+from dune.codegen.loopy.symbolic import FusedMultiplyAdd as FMA
 from dune.codegen.options import get_form_option
 from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
                                           EquidistantGeometryMixin,
@@ -15,12 +19,9 @@ from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
                                           restricted_name,
                                           name_cell_geometry
                                           )
-from dune.codegen.blockstructured.tools import (sub_element_inames,
-                                                name_point_in_macro,
-                                                )
+from dune.codegen.pdelab.tensors import name_matrix_inverse, name_determinant
+from dune.codegen.tools import get_pymbolic_basename
 from dune.codegen.ufl.modified_terminals import Restriction
-import pymbolic.primitives as prim
-from loopy.match import Writes
 
 
 @geometry_mixin("blockstructured_multilinear")
@@ -49,43 +50,25 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
                              prim.Power(get_form_option("number_of_blocks"), local_dimension()))
 
     def jacobian_determinant(self, o):
-        name = 'detjac'
-        self.define_jacobian_determinant(name)
-        return prim.Quotient(prim.Variable(name),
-                             prim.Power(get_form_option("number_of_blocks"), local_dimension()))
-
-    def define_jacobian_determinant(self, name):
-        temporary_variable(name, shape=(), managed=True)
-
-        determinant_signed = name_jacobian_determinant_signed(self)
+        jacobian = name_jacobian_matrix(self)
+        name = name_determinant(jacobian, (world_dimension(), world_dimension()), self)
 
-        return instruction(expression=prim.Call(prim.Variable("abs"), (prim.Variable(determinant_signed),)),
-                           assignee=prim.Variable(name),
-                           within_inames=frozenset(sub_element_inames() + self.quadrature_inames()),
-                           depends_on=frozenset({Writes(determinant_signed)})
-                           )
+        return prim.Quotient(prim.Call(prim.Variable("abs"), (prim.Variable(name),)),
+                             prim.Power(get_form_option("number_of_blocks"), local_dimension()))
 
     def jacobian_inverse(self, o):
-        restriction = enforce_boundary_restriction(self)
-
         assert(len(self.indices) == 2)
         i, j = self.indices
         self.indices = None
 
-        name = restricted_name("jit", restriction)
-        self.define_jacobian_inverse_transposed(name, restriction)
+        restriction = enforce_boundary_restriction(self)
+
+        jacobian = restricted_name(name_jacobian_matrix(self), restriction)
+        name = name_matrix_inverse(jacobian, (world_dimension(), world_dimension()), self)
 
         return prim.Product((prim.Subscript(prim.Variable(name), (j, i)),
                              get_form_option("number_of_blocks")))
 
-    def define_jacobian_inverse_transposed(self, name, restriction):
-        temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
-
-        jacobian = name_jacobian_matrix(self)
-        det_inv = name_jacobian_determinant_inverse(self)
-
-        compute_inverse_transposed(name, det_inv, jacobian, self)
-
 
 @geometry_mixin("blockstructured_axiparallel")
 class AxiparallelBlockStructuredGeometryMixin(AxiparallelGeometryMixin, BlockStructuredGeometryMixin):
@@ -167,12 +150,12 @@ def compute_jacobian(name, visitor):
         a, b, c = coefficients
 
         expr_jac = [None, None]
-        expr_jac[0] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (1,)),
-                                              prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)))),
-                                prim.Subscript(prim.Variable(b), (prim.Variable(jac_iname),))))
-        expr_jac[1] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (0,)),
-                                              prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)))),
-                                prim.Subscript(prim.Variable(c), (prim.Variable(jac_iname),))))
+        expr_jac[0] = FMA(prim.Subscript(pymbolic_pos, (1,)),
+                          prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)),
+                          prim.Subscript(prim.Variable(b), (prim.Variable(jac_iname),)))
+        expr_jac[1] = FMA(prim.Subscript(pymbolic_pos, (0,)),
+                          prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)),
+                          prim.Subscript(prim.Variable(c), (prim.Variable(jac_iname),)))
     elif world_dimension() == 3:
         a, b, c, d, e, f, g = coefficients
 
@@ -183,21 +166,20 @@ def compute_jacobian(name, visitor):
         # with k, l in {0,1,2} != i and k<l and vj = terms[i][j]
         for i in range(3):
             k, l = sorted(set(range(3)) - {i})
-            expr_jac[i] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (k,)), prim.Subscript(pymbolic_pos, (l,)),
-                                                  prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)))),
-                                    prim.Product((prim.Subscript(pymbolic_pos, (k,)),
-                                                  prim.Subscript(prim.Variable(terms[i][0]), (prim.Variable(jac_iname),)))),
-                                    prim.Product((prim.Subscript(pymbolic_pos, (l,)),
-                                                  prim.Subscript(prim.Variable(terms[i][1]), (prim.Variable(jac_iname),)))),
-                                    prim.Subscript(prim.Variable(terms[i][2]), (prim.Variable(jac_iname),))
-                                    ))
+            expr_jac[i] = FMA(prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)),
+                              prim.Subscript(pymbolic_pos, (k,)) * prim.Subscript(pymbolic_pos, (l,)),
+                              FMA(prim.Subscript(prim.Variable(terms[i][0]), (prim.Variable(jac_iname),)),
+                                  prim.Subscript(pymbolic_pos, (k,)),
+                                  FMA(prim.Subscript(prim.Variable(terms[i][1]), (prim.Variable(jac_iname),)),
+                                      prim.Subscript(pymbolic_pos, (l,)),
+                                      prim.Subscript(prim.Variable(terms[i][2]), (prim.Variable(jac_iname),)))))
     else:
         raise NotImplementedError()
 
     for i, expr in enumerate(expr_jac):
         instruction(expression=expr,
                     assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_iname), i)),
-                    within_inames=frozenset((jac_iname, ) + sub_element_inames() + visitor.quadrature_inames()),
+                    within_inames=frozenset((jac_iname, ) + visitor.quadrature_inames()),
                     depends_on=frozenset({Writes(get_pymbolic_basename(pymbolic_pos))} | {Writes(cd) for cd in coefficients})
                     )
 
@@ -213,121 +195,6 @@ def name_jacobian_matrix(visitor):
     return name
 
 
-def compute_determinant(name, matrix, visitor):
-    dim = world_dimension()
-    matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
-    if dim == 2:
-        expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
-                                     -1 * prim.Product((matrix_entry[1][0], matrix_entry[0][1]))))
-    elif dim == 3:
-        expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1], matrix_entry[2][2])),
-                                     prim.Product((matrix_entry[0][1], matrix_entry[1][2], matrix_entry[2][0])),
-                                     prim.Product((matrix_entry[0][2], matrix_entry[1][0], matrix_entry[2][1])),
-                                     -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1], matrix_entry[2][0])),
-                                     -1 * prim.Product((matrix_entry[0][0], matrix_entry[1][2], matrix_entry[2][1])),
-                                     -1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0], matrix_entry[2][2]))
-                                     ))
-    else:
-        raise NotImplementedError()
-    instruction(expression=expr_determinant,
-                assignee=prim.Variable(name),
-                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
-                depends_on=frozenset({Writes(matrix)})
-                )
-
-
-def define_jacobian_determinant(name, visitor):
-    temporary_variable(name, shape=(), managed=True)
-    jacobian = name_jacobian_matrix(visitor)
-    compute_determinant(name, jacobian, visitor)
-
-
-def define_jacobian_determinant_inverse(name, visitor):
-    temporary_variable(name, shape=(), managed=True)
-    determinant = name_jacobian_determinant_signed(visitor)
-    return instruction(expression=prim.Quotient(1., prim.Variable(determinant)),
-                       assignee=prim.Variable(name),
-                       within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
-                       depends_on=frozenset({Writes(determinant)})
-                       )
-
-
-def name_jacobian_determinant_signed(visitor):
-    name = "detjac_signed"
-    define_jacobian_determinant(name, visitor)
-    return name
-
-
-def name_jacobian_determinant_inverse(visitor):
-    name = "detjac_inverse"
-    define_jacobian_determinant_inverse(name, visitor)
-    return name
-
-
-def compute_inverse_transposed(name, det_inv, matrix, visitor):
-    dim = world_dimension()
-    matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
-    assignee = [[prim.Subscript(prim.Variable(name), (i, j)) for j in range(dim)] for i in range(dim)]
-    exprs = [[None for _ in range(dim)] for _ in range(dim)]
-
-    if dim == 2:
-        for i in range(2):
-            for j in range(2):
-                sign = 1. if i == j else -1.
-                exprs[i][j] = prim.Product((sign, prim.Variable(det_inv), matrix_entry[1 - i][1 - j]))
-    elif dim == 3:
-        exprs[0][0] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[1][1], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[1][2], matrix_entry[2][1]))))))
-        exprs[1][0] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[2][1]))))))
-        exprs[2][0] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[1][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1]))))))
-
-        exprs[0][1] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[1][2], matrix_entry[2][0]))))))
-        exprs[1][1] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[2][0]))))))
-        exprs[2][1] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][0]))))))
-
-        exprs[0][2] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][1])),
-                                              -1 * prim.Product((matrix_entry[1][1], matrix_entry[2][0]))))))
-        exprs[1][2] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][1])),
-                                              -1 * prim.Product((matrix_entry[0][1], matrix_entry[2][0]))))))
-        exprs[2][2] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
-                                              -1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0]))))))
-    else:
-        raise NotImplementedError
-    for j in range(dim):
-        for i in range(dim):
-            instruction(expression=exprs[i][j],
-                        assignee=assignee[i][j],
-                        within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
-                        depends_on=frozenset({Writes(matrix), Writes(det_inv)}))
-
-
-def define_jacobian_inverse_transposed(name, visitor):
-    temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
-    jacobian = name_jacobian_matrix(visitor)
-    det_inv = name_jacobian_determinant_inverse(visitor)
-    compute_inverse_transposed(name, det_inv, jacobian, visitor)
-
-
-def name_jacobian_inverse_transposed(restriction, visitor):
-    name = restricted_name("jit", restriction)
-    define_jacobian_inverse_transposed(name, visitor)
-    return name
-
-
 def compute_multilinear_to_global_transformation(name, local, visitor):
     dim = world_dimension()
     temporary_variable(name, shape=(dim,), managed=True)
@@ -347,19 +214,22 @@ def compute_multilinear_to_global_transformation(name, local, visitor):
     # global[d] = T(local)[d]
     if dim == 2:
         a_pym, b_pym, c_pym = coeffs_pym
-        expr = a_pym * local_pym[0] * local_pym[1] + b_pym * local_pym[0] + c_pym * local_pym[1] + corner_0_pym
+        expr = FMA(a_pym, local_pym[0] * local_pym[1], FMA(b_pym, local_pym[0], FMA(c_pym, local_pym[1], corner_0_pym)))
     elif dim == 3:
         a_pym, b_pym, c_pym, d_pym, e_pym, f_pym, g_pym = coeffs_pym
-        expr = (a_pym * local_pym[0] * local_pym[1] * local_pym[2] + b_pym * local_pym[0] * local_pym[1] +
-                c_pym * local_pym[0] * local_pym[2] + d_pym * local_pym[1] * local_pym[2] +
-                e_pym * local_pym[0] + f_pym * local_pym[1] + g_pym * local_pym[2] + corner_0_pym)
+        expr = FMA(a_pym * local_pym[0], local_pym[1] * local_pym[2],
+                   FMA(b_pym, local_pym[0] * local_pym[1],
+                       FMA(c_pym, local_pym[0] * local_pym[2],
+                           FMA(d_pym, local_pym[1] * local_pym[2],
+                               FMA(e_pym, local_pym[0],
+                                   FMA(f_pym, local_pym[1], FMA(g_pym, local_pym[2], corner_0_pym)))))))
     else:
         raise NotImplementedError
 
     assignee = prim.Subscript(prim.Variable(name), (dim_pym,))
 
     instruction(assignee=assignee, expression=expr,
-                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames() + (dim_pym.name,)),
+                within_inames=frozenset(visitor.quadrature_inames() + (dim_pym.name,)),
                 within_inames_is_final=True,
                 depends_on=frozenset({Writes(get_pymbolic_basename(local)), Writes(corners)})
                 )
@@ -376,13 +246,14 @@ def compute_axiparallel_to_global_transformation(name, local, visitor):
     dim_pym = prim.Variable(component_iname('to_global'))
 
     # global[d] = lower_left[d] + local[d] * (upper_right[d] - lower_left[d])
-    expr = (prim.Subscript(prim.Variable(corners), (0, dim_pym)) +
-            prim.Subscript(local, (dim_pym,)) * (prim.Subscript(prim.Variable(corners), (2**dim - 1, dim_pym)) -
-                                                 prim.Subscript(prim.Variable(corners), (0, dim_pym))))
+    expr = FMA(prim.Subscript(prim.Variable(corners), (2**dim - 1, dim_pym)) -
+               prim.Subscript(prim.Variable(corners), (0, dim_pym)),
+               prim.Subscript(local, (dim_pym,)), prim.Subscript(prim.Variable(corners), (0, dim_pym)))
+
     assignee = prim.Subscript(prim.Variable(name), (dim_pym,))
 
     instruction(assignee=assignee, expression=expr,
-                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames() + (dim_pym.name,)),
+                within_inames=frozenset(visitor.quadrature_inames() + (dim_pym.name,)),
                 within_inames_is_final=True,
                 depends_on=frozenset({Writes(get_pymbolic_basename(local)), Writes(corners)})
                 )

python/dune/codegen/blockstructured/quadrature.py

 from dune.codegen.error import CodegenError
 from dune.codegen.generation import quadrature_mixin
 from dune.codegen.pdelab.quadrature import GenericQuadratureMixin
-from dune.codegen.blockstructured.tools import name_point_in_macro
+from dune.codegen.blockstructured.tools import name_point_in_macro, sub_element_inames, remove_sub_element_inames
 import pymbolic.primitives as prim
 
 
 @quadrature_mixin("blockstructured")
 class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
+    @staticmethod
+    def _subscript_without_sub_element_inames(s):
+        return prim.Subscript(s.aggregate, remove_sub_element_inames(s.index_tuple))
+
     def quadrature_position_in_micro(self, index=None):
-        return GenericQuadratureMixin.quadrature_position(self, index)
+        return self._subscript_without_sub_element_inames(super().quadrature_position(index))
 
     def quadrature_position_in_macro(self, index=None):
         original = self.quadrature_position_in_micro(index)
@@ -20,3 +24,9 @@ class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
 
     def quadrature_position(self, index=None):
         raise CodegenError('Decide if the quadrature point is in the macro or micro element')
+
+    def quadrature_inames(self):
+        return super().quadrature_inames() + sub_element_inames()
+
+    def quadrature_weight(self, o):
+        return self._subscript_without_sub_element_inames(super().quadrature_weight(o))

python/dune/codegen/blockstructured/tools.py

@@ -22,6 +22,11 @@ def sub_element_inames():
     return inames
 
 
+def remove_sub_element_inames(indices):
+    assert isinstance(indices, tuple)
+    return tuple(set(indices) - set(sub_element_inames()) - set(prim.Variable(i) for i in sub_element_inames()))
+
+
 # compute sequential index for given tensor index, the same as index in base-k to base-10
 def tensor_index_to_sequential_index(indices, k):
     return prim.Sum(tuple(prim.Variable(index) * k ** i for i, index in enumerate(indices)))
@@ -64,7 +69,7 @@ def define_point_in_macro(name, point_in_micro, visitor):
         # TODO relax within inames
         instruction(assignee=prim.Subscript(prim.Variable(name), (i,)),
                     expression=expr,
-                    within_inames=frozenset(subelem_inames + visitor.quadrature_inames()),
+                    within_inames=frozenset(visitor.quadrature_inames()),
                     tags=frozenset({subelem_inames[i]})
                     )
 

python/dune/codegen/pdelab/tensors.py

 """ Code generation for explicitly specified tensors """
 
 from dune.codegen.generation import (get_counted_variable,
-                                     domain,
                                      kernel_cached,
-                                     iname,
                                      instruction,
                                      temporary_variable,
                                      )
+from dune.codegen.loopy.symbolic import FusedMultiplyAdd as FMA
+from loopy.match import Writes
 
 import pymbolic.primitives as prim
 import numpy as np
@@ -14,6 +14,116 @@ import loopy as lp
 import itertools as it
 
 
+def define_determinant(name, matrix, shape, visitor):
+    temporary_variable(name, managed=True)
+
+    assert len(shape) == 2 and shape[0] == shape[1]
+    dim = shape[0]
+
+    matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
+    if dim == 2:
+        expr_determinant = FMA(matrix_entry[0][0], matrix_entry[1][1], -1 * matrix_entry[1][0] * matrix_entry[0][1])
+
+    elif dim == 3:
+        fma_A = FMA(matrix_entry[1][1], matrix_entry[2][2], -1 * matrix_entry[1][2] * matrix_entry[2][1])
+        fma_B = FMA(matrix_entry[1][0], matrix_entry[2][2], -1 * matrix_entry[1][2] * matrix_entry[2][0])
+        fma_C = FMA(matrix_entry[1][0], matrix_entry[2][1], -1 * matrix_entry[1][1] * matrix_entry[2][0])
+
+        expr_determinant = FMA(matrix_entry[0][2], fma_C,
+                               FMA(matrix_entry[0][0], fma_A, -1 * matrix_entry[0][1] * fma_B))
+    else:
+        raise NotImplementedError()
+    instruction(expression=expr_determinant,
+                assignee=prim.Variable(name),
+                within_inames=frozenset(visitor.quadrature_inames()),
+                depends_on=frozenset({Writes(matrix)})
+                )
+
+
+def define_determinant_inverse(name, matrix, shape, visitor):
+    det = name_determinant(matrix, shape, visitor)
+
+    temporary_variable(name, managed=True)
+
+    instruction(expression=prim.Quotient(1, prim.Variable(det)),
+                assignee=prim.Variable(name),
+                within_inames=frozenset(visitor.quadrature_inames()),
+                depends_on=frozenset({Writes(matrix), Writes(det)})
+                )
+
+
+def define_matrix_inverse(name, name_inv, shape, visitor):
+    temporary_variable(name_inv, shape=shape, managed=True)
+
+    det_inv = name_determinant_inverse(name, shape, visitor)
+
+    assert len(shape) == 2 and shape[0] == shape[1]
+    dim = shape[0]
+
+    matrix_entry = [[prim.Subscript(prim.Variable(name), (i, j)) for j in range(dim)] for i in range(dim)]
+    assignee = [[prim.Subscript(prim.Variable(name_inv), (i, j)) for j in range(dim)] for i in range(dim)]
+    exprs = [[None for _ in range(dim)] for _ in range(dim)]
+
+    if dim == 2:
+        for i in range(2):
+            for j in range(2):
+                sign = 1. if i == j else -1.
+                exprs[i][j] = prim.Product((sign, prim.Variable(det_inv), matrix_entry[1 - i][1 - j]))
+    elif dim == 3:
+        exprs[0][0] = prim.Variable(det_inv) * FMA(matrix_entry[1][1], matrix_entry[2][2],
+                                                   -1 * matrix_entry[1][2] * matrix_entry[2][1])
+        exprs[1][0] = prim.Variable(det_inv) * FMA(matrix_entry[0][1], matrix_entry[2][2],
+                                                   -1 * matrix_entry[0][2] * matrix_entry[2][1]) * -1
+        exprs[2][0] = prim.Variable(det_inv) * FMA(matrix_entry[0][1], matrix_entry[1][2],
+                                                   -1 * matrix_entry[0][2] * matrix_entry[1][1])
+
+        exprs[0][1] = prim.Variable(det_inv) * FMA(matrix_entry[1][0], matrix_entry[2][2],
+                                                   -1 * matrix_entry[1][2] * matrix_entry[2][0]) * -1
+        exprs[1][1] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[2][2],
+                                                   -1 * matrix_entry[0][2] * matrix_entry[2][0])
+        exprs[2][1] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[1][2],
+                                                   -1 * matrix_entry[0][2] * matrix_entry[1][0]) * -1
+
+        exprs[0][2] = prim.Variable(det_inv) * FMA(matrix_entry[1][0], matrix_entry[2][1],
+                                                   -1 * matrix_entry[1][1] * matrix_entry[2][0])
+        exprs[1][2] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[2][1],
+                                                   -1 * matrix_entry[0][1] * matrix_entry[2][0]) * -1
+        exprs[2][2] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[1][1],
+                                                   -1 * matrix_entry[0][1] * matrix_entry[1][0])
+    else:
+        raise NotImplementedError
+    for j in range(dim):
+        for i in range(dim):
+            instruction(expression=exprs[i][j],
+                        assignee=assignee[i][j],
+                        within_inames=frozenset(visitor.quadrature_inames()),
+                        depends_on=frozenset({Writes(name), Writes(det_inv)}))
+
+
+def name_determinant(matrix, shape, visitor):
+    name = matrix + "_det"
+
+    define_determinant(name, matrix, shape, visitor)
+
+    return name
+
+
+def name_determinant_inverse(matrix, shape, visitor):
+    name = matrix + "_det_inv"
+
+    define_determinant_inverse(name, matrix, shape, visitor)
+
+    return name
+
+
+def name_matrix_inverse(name, shape, visitor):
+    name_inv = name + "_inv"
+
+    define_matrix_inverse(name, name_inv, shape, visitor)
+
+    return name_inv
+
+
 def define_assembled_tensor(name, expr, visitor):
     temporary_variable(name,
                        shape=expr.ufl_shape,
@@ -64,16 +174,19 @@ def pymbolic_matrix_inverse(o, visitor):
 
     # If code generation time inversion failed, we assemble it in C++
     # and invert it there.
-    name = name_assembled_tensor(o.ufl_operands[0], visitor)
-
-    instruction(code="{}.invert();".format(name),
-                within_inames=frozenset(visitor.quadrature_inames()),
-                depends_on=frozenset({lp.match.Writes(name),
-                                      lp.match.Tagged("sumfact_stage1"),
-                                      }),
-                tags=frozenset({"quad"}),
-                assignees=frozenset({name}),
-                )
+    expr = o.ufl_operands[0]
+    name = name_assembled_tensor(expr, visitor)
+
+    if expr.shape[0] <= 3:
+        name = name_matrix_inverse(name, expr.ufl_shape, visitor)
+    else:
+        instruction(code="{}.invert();".format(name),
+                    within_inames=frozenset(visitor.quadrature_inames()),
+                    depends_on=frozenset({lp.match.Writes(name),
+                                          lp.match.Tagged("sumfact_stage1"),
+                                          }),
+                    tags=frozenset({name}),
+                    )
 
     visitor.indices = indices
     return prim.Variable(name)