[!263] Implement code generation for matrix inversion in C++

Merge branch 'feature/matrix-inversion' into 'master'

ref:dominic/dune-perftool Matrix inversion at code generation time does only
work to a very limited extent (up to n=4). We can instead assemble the tensor
in C++ and invert it there (e.g. using Dune::FieldMatrix)

This fixes [#123].

Still TODO:

-   \[ \] Vectorized Inversion

See merge request [dominic/dune-perftool!263]

  [#123]: gitlab.dune-project.org/NoneNone/issues/123
  [dominic/dune-perftool!263]: gitlab.dune-project.org/dominic/dune-perftool/merge_requests/263


Closes #123

python/dune/perftool/pdelab/__init__.py

@@ -32,7 +32,10 @@ from dune.perftool.pdelab.quadrature import (pymbolic_quadrature_weight,
                                              )
 from dune.perftool.pdelab.spaces import (lfs_inames,
                                          )
-from dune.perftool.pdelab.tensors import pymbolic_list_tensor, pymbolic_identity
+from dune.perftool.pdelab.tensors import (pymbolic_list_tensor,
+                                          pymbolic_identity,
+                                          pymbolic_matrix_inverse,
+                                          )
 
 
 class PDELabInterface(object):
@@ -112,6 +115,9 @@ class PDELabInterface(object):
     def pymbolic_identity(self, o):
         return pymbolic_identity(o)
 
+    def pymbolic_matrix_inverse(self, o):
+        return pymbolic_matrix_inverse(o, self.visitor)
+
     #
     # Geometry related generator functions
     #

python/dune/perftool/pdelab/tensors.py

@@ -11,6 +11,7 @@ from dune.perftool.generation import (get_counted_variable,
 import pymbolic.primitives as prim
 import numpy as np
 import loopy as lp
+import itertools as it
 
 
 def define_list_tensor(name, expr, visitor, stack=()):
@@ -66,3 +67,42 @@ def pymbolic_identity(expr):
                        )
     define_identity(name, expr)
     return prim.Variable(name)
+
+
+def define_assembled_tensor(name, expr, visitor):
+    temporary_variable(name,
+                       shape=expr.ufl_shape,
+                       shape_impl=('fm',))
+    for indices in it.product(*tuple(range(i) for i in expr.ufl_shape)):
+        visitor.indices = indices
+        instruction(assignee=prim.Subscript(prim.Variable(name), indices),
+                    expression=visitor.call(expr),
+                    forced_iname_deps=frozenset(visitor.interface.quadrature_inames()),
+                    depends_on=frozenset({lp.match.Tagged("sumfact_stage1")}),
+                    tags=frozenset({"quad"}),
+                    )
+
+
+@kernel_cached
+def name_assembled_tensor(o, visitor):
+    name = get_counted_variable("assembled_tensor")
+    define_assembled_tensor(name, o, visitor)
+    return name
+
+
+@kernel_cached
+def pymbolic_matrix_inverse(o, visitor):
+    indices = visitor.indices
+    visitor.indices = None
+    name = name_assembled_tensor(o.ufl_operands[0], visitor)
+
+    instruction(code="{}.invert();".format(name),
+                within_inames=frozenset(visitor.interface.quadrature_inames()),
+                depends_on=frozenset({lp.match.Writes(name),
+                                      lp.match.Tagged("sumfact_stage1"),
+                                      }),
+                tags=frozenset({"quad"}),
+                )
+
+    visitor.indices = indices
+    return prim.Variable(name)

python/dune/perftool/ufl/preprocess.py

@@ -3,6 +3,10 @@
 from dune.perftool.generation import run_hook, ReturnArg
 import ufl.classes as uc
 import ufl.algorithms.apply_function_pullbacks as afp
+import ufl.algorithms.apply_algebra_lowering as aal
+import ufl.algorithms.apply_derivatives as ad
+
+from dune.perftool.options import get_form_option
 
 from pytools import memoize
 
@@ -35,7 +39,8 @@ def preprocess_form(form):
                                                           uc.FacetNormal,
                                                           uc.JacobianDeterminant,
                                                           uc.JacobianInverse,
-                                                          )
+                                                          ),
+                                 do_estimate_degrees=not get_form_option("quadrature_order"),
                                  )
 
     formdata.preprocessed_form = apply_default_transformations(formdata.preprocessed_form)
@@ -58,3 +63,11 @@ def apply_default_transformations(form):
     form = transform_form(form, pushdown_indexed)
 
     return form
+
+
+# Monkey patch UFL, such that we invert matrices in C++ instead of Python.
+# The latter only works for very small matrices. If this causes a problem at
+# some point, we should guard this monkey patch with an option.
+
+aal.LowerCompoundAlgebra.inverse = lambda s, o, A: s.reuse_if_untouched(o, A)
+ad.GenericDerivativeRuleset.inverse = lambda s, o, A: -uc.Dot(uc.Dot(o, A), o)

python/dune/perftool/ufl/visitor.py

@@ -267,6 +267,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
     def identity(self, o):
         return self.interface.pymbolic_identity(o)
 
+    def inverse(self, o):
+        return self.interface.pymbolic_matrix_inverse(o)
+
     #
     # Handlers for arithmetic operators and functions
     # Those handlers would be valid in any code going from UFL to pymbolic