Point evaluation fails with 64bit integers.

[dham]

If Firedrake is installed with --petsc-int-type int64 then point evaluation generates illegal code which causes a compiler error. MFE:

from firedrake import *

mesh = UnitSquareMesh(2, 2)

fs = FunctionSpace(mesh, "CG", 1)
f = Function(fs)

f.at([0,0])

The issue appears to be that wrap_to_reference_coords takes int64_t arguments but is passed arguments of type long.

[dham]

I have traced this to an issue with how loopy maps numpy types to c types. I think I have a solution for that, but it requires an extra header file.

[wence-]

Can we just fix it upstream?

[kaushikcfd]

@dham: One way I could think of including the header is --

index 25b19080..9e05446a 100644
--- a/loopy/target/c/__init__.py
+++ b/loopy/target/c/__init__.py
@@ -77,6 +77,11 @@ class DTypeRegistryWrapper(object):

 # {{{ preamble generator

+def c99_preamble_generator(preamble_info):
+    if any(dtype.is_integral() for dtype in preamble_info.seen_dtypes):
+        yield("10_inttypes", "#include <inttypes.h>")
+
+
 def _preamble_generator(preamble_info):
     integer_type_names = ["int8", "int16", "int32", "int64"]

@@ -1080,4 +1085,21 @@ def generate_header(kernel, codegen_result=None):

 # }}}

+
+# {{{
+
+class C99Target(CTarget):
+    def get_device_ast_builder(self):
+        return C99ASTBuilder(self)
+
+
+class C99ASTBuilder(CASTBuilder):
+    def preamble_generators(self):
+        return (
+                super(CASTBuilder, self).preamble_generators() + [
+                    c99_preamble_generator,
+                    ])
+
+# }}}

Directly including the preamble in CASTBuilder would be a bit problematic as OpenCLCASTBuilder, CUDACASTBuilder, etc. derive from CASTBuilder. (cc @inducer)

[inducer]

I broadly agree with @kaushikcfd. One possibility would be to reshuffle the inheritance structure of the C targets so as to distinguish the ones from which we expect full C99 compliance.

[kaushikcfd]

Related MR: https://gitlab.tiker.net/inducer/loopy/merge_requests/357/diffs

[dham]

I think I don't understand that MR. Is there an unpushed change to compete that I should be looking for?

[kaushikcfd]

There is also a change in compyte: https://github.com/inducer/compyte/pull/28/files

[dham]

OK, That makes perfect sense now. If those land then we can update the Firedrake loopy branch accordingly.

[inducer]

They're on the way, just pending a final look by Kaushik. (Thanks for working on this, @kaushikcfd!)

[danshapero]

I've just built Firedrake with 64-bit indices and it looks like this is still a problem. But the Firedrake fork of loopy hasn't yet included the changes that @kaushikcfd and @inducer added in October. I tried to merge in upstream loopy but got loads of conflicts.

[kaushikcfd]

But the Firedrake fork of loopy hasn't yet included the changes that @kaushikcfd and @inducer added in October

These changes were merged into firedrakeproject/loopy. See: firedrakeproject/loopy@c91cb1b

[wence-]

@danshapero when you have time can you confirm things now work (I think they should).

[danshapero]

I just installed again with 64-bit indices and ran David's failing example and I'm still getting a compilation error on my system. In the generated code, it looks like there's a forward declaration of a function wrap_evaluate that conflicts with the actual definition. I've attached a zip archive with the generated code, the logfile, and the compile errors.
generated-code.zip

[sv2518]

Maybe this would be the right thing to do?

diff --git a/firedrake/pointeval_utils.py b/firedrake/pointeval_utils.py
index 07923496..ade22938 100644
--- a/firedrake/pointeval_utils.py
+++ b/firedrake/pointeval_utils.py
@@ -125,7 +125,7 @@ def compile_element(expression, coordinates, parameters=None):
         code["map_args"] = "f->coords_map, f->f_map"
 
     evaluate_template_c = """
-static inline void wrap_evaluate(%(scalar_type)s* const result, %(scalar_type)s* const X, int const start, int const end%(layers_arg)s,
+static inline void wrap_evaluate(%(scalar_type)s* const result, %(scalar_type)s* const X, %(IntType)s const start, %(IntType)s const end%(layers_arg)s,
     %(scalar_type)s const *__restrict__ coords, %(scalar_type)s const *__restrict__ f, %(wrapper_map_args)s);

 int evaluate(struct Function *f, %(scalar_type)s *x, %(scalar_type)s *result)

[danshapero]

That worked perfectly, thanks @sv2518 !