diff --git a/astnodes.py b/astnodes.py
index 5a3602f410570a1978d946c2671342e3f25d91dd..f6c2d24de3cf93660b49d1ca563e3e92d49b694e 100644
--- a/astnodes.py
+++ b/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,
diff --git a/cpu/cpujit.py b/cpu/cpujit.py
index 18429d1200c9a04f0c50cb7abc994c43398ff63a..2bc08985668322a96bba9a4a737d9d98bb4ccc57 100644
--- a/cpu/cpujit.py
+++ b/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
diff --git a/field.py b/field.py
index 18ea09dafc997cec4f129fb507e278bd0e0f8465..2b3c19e9da4a77dd9d718bb31d9d8f6395ceb1cd 100644
--- a/field.py
+++ b/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):
diff --git a/transformations.py b/transformations.py
index c93d98ca291521838efda267d9e62bbe4a684ecf..1d2dfd5244a4297ed56b9f40e0f2211c0e207402 100644
--- a/transformations.py
+++ b/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)