pystencils: additional checks when calling kernel

- check that fixed size kernels are called with arrays of the correct size
- checks that layout of compilation matches runtime layout
- not allowed any more to mix fixed & and variable sized fields in a kernel

astnodes.py

@@ -40,9 +40,10 @@ class Node(object):
 class KernelFunction(Node):
 
     class Argument:
-        def __init__(self, name, dtype):
+        def __init__(self, name, dtype, kernelFunctionNode):
+            from pystencils.transformations import symbolNameToVariableName
             self.name = name
-            self.dtype = dtype  # TODO ordentliche Klasse
+            self.dtype = dtype
             self.isFieldPtrArgument = False
             self.isFieldShapeArgument = False
             self.isFieldStrideArgument = False
@@ -63,6 +64,11 @@ class KernelFunction(Node):
                 self.isFieldArgument = True
                 self.fieldName = name[len(Field.STRIDE_PREFIX):]
 
+            self.field = None
+            if self.isFieldArgument:
+                fieldMap = {symbolNameToVariableName(f.name): f for f in kernelFunctionNode.fieldsAccessed}
+                self.field = fieldMap[self.fieldName]
+
         def __repr__(self):
             return '<{0} {1}>'.format(self.dtype, self.name)
 
@@ -105,12 +111,11 @@ class KernelFunction(Node):
 
     def _updateParameters(self):
         undefinedSymbols = self._body.undefinedSymbols - self.globalVariables
-        self._parameters = [KernelFunction.Argument(s.name, s.dtype) for s in undefinedSymbols]
+        self._parameters = [KernelFunction.Argument(s.name, s.dtype, self) for s in undefinedSymbols]
         self._parameters.sort(key=lambda l: (l.fieldName, l.isFieldPtrArgument, l.isFieldShapeArgument,
                                              l.isFieldStrideArgument, l.name),
                               reverse=True)
 
-
     def __str__(self):
         self._updateParameters()
         return '{0} {1}({2})\n{3}'.format(type(self).__name__, self.functionName, self.parameters,

cpu/cpujit.py

@@ -10,7 +10,7 @@ import hashlib
 from pystencils.transformations import symbolNameToVariableName
 
 CONFIG_GCC = {
-    'compiler': 'g++',
+    'compiler': 'g++-4.8',
     'flags': '-Ofast -DNDEBUG -fPIC -shared -march=native -fopenmp',
 }
 CONFIG_INTEL = {
@@ -108,11 +108,24 @@ def compileAndLoad(kernelFunctionNode):
 def buildCTypeArgumentList(parameterSpecification, argumentDict):
     argumentDict = {symbolNameToVariableName(k): v for k, v in argumentDict.items()}
     ctArguments = []
+    arrayShapes = set()
     for arg in parameterSpecification:
         if arg.isFieldArgument:
             field = argumentDict[arg.fieldName]
+            symbolicField = arg.field
             if arg.isFieldPtrArgument:
                 ctArguments.append(field.ctypes.data_as(ctypeFromString(arg.dtype)))
+                if symbolicField.hasFixedShape:
+                    if tuple(int(i) for i in symbolicField.shape) != field.shape:
+                        raise ValueError("Passed array '%s' has shape %s which does not match expected shape %s" %
+                                         (arg.fieldName, str(field.shape), str(symbolicField.shape)))
+                if symbolicField.hasFixedShape:
+                    if tuple(int(i) * field.itemsize for i in symbolicField.strides) != field.strides:
+                        raise ValueError("Passed array '%s' has strides %s which does not match expected strides %s" %
+                                         (arg.fieldName, str(field.strides), str(symbolicField.strides)))
+
+                if not symbolicField.isIndexField:
+                    arrayShapes.add(field.shape[:symbolicField.spatialDimensions])
             elif arg.isFieldShapeArgument:
                 dataType = ctypeFromString(arg.dtype, includePointers=False)
                 ctArguments.append(field.ctypes.shape_as(dataType))
@@ -130,6 +143,9 @@ def buildCTypeArgumentList(parameterSpecification, argumentDict):
             param = argumentDict[arg.name]
             expectedType = ctypeFromString(arg.dtype)
             ctArguments.append(expectedType(param))
+
+    if len(arrayShapes) > 1:
+        raise ValueError("All passed arrays have to have the same size " + str(arrayShapes))
     return ctArguments
 
 

field.py

@@ -46,8 +46,34 @@ class Field:
         Eq(dst_C^0, src_C^0)
         Eq(dst_C^1, src_S^1)
         Eq(dst_C^2, src_N^2)
-
     """
+
+    @staticmethod
+    def createGeneric(fieldName, spatialDimensions, dtype=np.float64, indexDimensions=0, layout='numpy'):
+        """
+        Creates a generic field where the field size is not fixed i.e. can be called with arrays of different sizes
+
+        :param fieldName: symbolic name for the field
+        :param dtype: numpy data type of the array the kernel is called with later
+        :param spatialDimensions: see documentation of Field
+        :param indexDimensions: see documentation of Field
+        :param layout: tuple specifying the loop ordering of the spatial dimensions e.g. (2, 1, 0 ) means that
+                       the outer loop loops over dimension 2, the second outer over dimension 1, and the inner loop
+                       over dimension 0. Also allowed: the strings 'numpy' (0,1,..d) or 'reverseNumpy' (d, ..., 1, 0)
+        """
+        if isinstance(layout, str) and (layout == 'numpy' or layout.lower() == 'c'):
+            layout = tuple(range(spatialDimensions))
+        elif isinstance(layout, str) and (layout == 'reverseNumpy' or layout.lower() == 'f'):
+            layout = tuple(reversed(range(spatialDimensions)))
+        if len(layout) != spatialDimensions:
+            raise ValueError("Layout")
+        shapeSymbol = IndexedBase(TypedSymbol(Field.SHAPE_PREFIX + fieldName, Field.SHAPE_DTYPE), shape=(1,))
+        strideSymbol = IndexedBase(TypedSymbol(Field.STRIDE_PREFIX + fieldName, Field.STRIDE_DTYPE), shape=(1,))
+        totalDimensions = spatialDimensions + indexDimensions
+        shape = tuple([shapeSymbol[i] for i in range(totalDimensions)])
+        strides = tuple([strideSymbol[i] for i in range(totalDimensions)])
+        return Field(fieldName, dtype, layout, shape, strides)
+
     @staticmethod
     def createFromNumpyArray(fieldName, npArray, indexDimensions=0):
         """
@@ -66,43 +92,43 @@ class Field:
         assert len(spatialLayout) == spatialDimensions
 
         strides = tuple([s // np.dtype(npArray.dtype).itemsize for s in npArray.strides])
-        shape = tuple([int(s) for s in npArray.shape])
+        shape = tuple(int(s) for s in npArray.shape)
 
         return Field(fieldName, npArray.dtype, spatialLayout, shape, strides)
 
     @staticmethod
-    def createGeneric(fieldName, spatialDimensions, dtype=np.float64, indexDimensions=0, layout='numpy'):
+    def createFixedSize(fieldName, shape, indexDimensions=0, dtype=np.float64, layout='numpy'):
         """
-        Creates a generic field where the field size is not fixed i.e. can be called with arrays of different sizes
+        Creates a field with fixed sizes i.e. can be called only wity arrays of the same size and layout
 
         :param fieldName: symbolic name for the field
+        :param shape: overall shape of the array
+        :param indexDimensions: how many of the trailing dimensions are interpreted as index (as opposed to spatial)
         :param dtype: numpy data type of the array the kernel is called with later
-        :param spatialDimensions: see documentation of Field
-        :param indexDimensions: see documentation of Field
-        :param layout: tuple specifying the loop ordering of the spatial dimensions e.g. (2, 1, 0 ) means that
-                       the outer loop loops over dimension 2, the second outer over dimension 1, and the inner loop
-                       over dimension 0. Also allowed: the strings 'numpy' (0,1,..d) or 'reverseNumpy' (d, ..., 1, 0)
+        :param layout: see createGeneric
         """
+        spatialDimensions = len(shape) - indexDimensions
+        assert spatialDimensions >= 1
+
         if isinstance(layout, str) and (layout == 'numpy' or layout.lower() == 'c'):
             layout = tuple(range(spatialDimensions))
         elif isinstance(layout, str) and (layout == 'reverseNumpy' or layout.lower() == 'f'):
             layout = tuple(reversed(range(spatialDimensions)))
-        if len(layout) != spatialDimensions:
-            raise ValueError("Layout")
-        shapeSymbol = IndexedBase(TypedSymbol(Field.SHAPE_PREFIX + fieldName, Field.SHAPE_DTYPE), shape=(1,))
-        strideSymbol = IndexedBase(TypedSymbol(Field.STRIDE_PREFIX + fieldName, Field.STRIDE_DTYPE), shape=(1,))
-        totalDimensions = spatialDimensions + indexDimensions
-        shape = tuple([shapeSymbol[i] for i in range(totalDimensions)])
-        strides = tuple([strideSymbol[i] for i in range(totalDimensions)])
-        return Field(fieldName, dtype, layout, shape, strides)
+
+        shape = tuple(int(s) for s in shape)
+        strides = computeStrides(shape, layout)
+        return Field(fieldName, dtype, layout[:spatialDimensions], shape, strides)
 
     def __init__(self, fieldName, dtype, layout, shape, strides):
         """Do not use directly. Use static create* methods"""
         self._fieldName = fieldName
         self._dtype = numpyDataTypeToC(dtype)
-        self._layout = layout
+        self._layout = normalizeLayout(layout)
         self.shape = shape
         self.strides = strides
+        # index fields are currently only used for boundary handling
+        # the coordinates are not the loop counters in that case, but are read from this index field
+        self.isIndexField = False
 
     @property
     def spatialDimensions(self):
@@ -124,6 +150,14 @@ class Field:
     def spatialShape(self):
         return self.shape[:self.spatialDimensions]
 
+    @property
+    def hasFixedShape(self):
+        try:
+            [int(i) for i in self.shape]
+            return True
+        except TypeError:
+            return False
+
     @property
     def indexShape(self):
         return self.shape[self.spatialDimensions:]
@@ -308,7 +342,34 @@ def getLayoutFromNumpyArray(arr, indexDimensionIds=[]):
     """
     coordinates = list(range(len(arr.shape)))
     relevantStrides = [stride for i, stride in enumerate(arr.strides) if i not in indexDimensionIds]
-    return tuple(x for (y, x) in sorted(zip(relevantStrides, coordinates), key=lambda pair: pair[0], reverse=True))
+    result = [x for (y, x) in sorted(zip(relevantStrides, coordinates), key=lambda pair: pair[0], reverse=True)]
+    return normalizeLayout(result)
+
+
+def normalizeLayout(layout):
+    """Takes a layout tuple and subtracts the minimum from all entries"""
+    minEntry = min(layout)
+    return tuple(i - minEntry for i in layout)
+
+
+def computeStrides(shape, layout):
+    """
+    Computes strides assuming no padding exists
+    :param shape: shape (size) of array
+    :param layout: layout specification as tuple
+    :return: strides in elements, not in bytes
+    """
+    layout = list(reversed(layout))
+    N = len(shape)
+    assert len(layout) == N
+    assert len(set(layout)) == N
+    strides = [0] * N
+    product = 1
+    for i in range(N):
+        j = layout.index(i)
+        strides[j] = product
+        product *= shape[j]
+    return tuple(strides)
 
 
 def numpyDataTypeToC(dtype):

transformations.py

@@ -43,15 +43,24 @@ def makeLoopOverDomain(body, functionName, iterationSlice=None, ghostLayers=None
     if loopOrder is None:
         loopOrder = getOptimalLoopOrdering(fields)
 
-    shapes = set([f.spatialShape for f in fields])
-
-    if len(shapes) > 1:
-        nrOfFixedSizedFields = 0
-        for shape in shapes:
-            if not isinstance(shape[0], sp.Basic):
-                nrOfFixedSizedFields += 1
-        assert nrOfFixedSizedFields <= 1, "Differently sized field accesses in loop body: " + str(shapes)
-    shape = list(shapes)[0]
+    nrOfFixedShapedFields = 0
+    for f in fields:
+        if f.hasFixedShape:
+            nrOfFixedShapedFields += 1
+
+    if nrOfFixedShapedFields > 0 and nrOfFixedShapedFields != len(fields):
+        fixedFieldNames = ",".join([f.name for f in fields if f.hasFixedShape])
+        varFieldNames = ",".join([f.name for f in fields if not f.hasFixedShape])
+        msg = "Mixing fixed-shaped and variable-shape fields in a single kernel is not possible\n"
+        msg += "Variable shaped: %s \nFixed shaped:    %s" % (varFieldNames, fixedFieldNames)
+        raise ValueError(msg)
+
+    shapeSet = set([f.spatialShape for f in fields])
+    if nrOfFixedShapedFields == len(fields):
+        if len(shapeSet) != 1:
+            raise ValueError("Differently sized field accesses in loop body: " + str(shapeSet))
+
+    shape = list(shapeSet)[0]
 
     if iterationSlice is not None:
         iterationSlice = normalizeSlice(iterationSlice, shape)