Commit 39a096fa authored by Jan Hoenig's avatar Jan Hoenig
Browse files

merged

parents 2d654ff4 c8b455fe
import sympy as sp
from sympy.tensor import IndexedBase
from pystencils.field import Field
from pystencils.types import TypedSymbol, DataType, get_type_from_sympy, _c_dtype_dict
from pystencils.types import TypedSymbol, createType, get_type_from_sympy
class Node(object):
......@@ -261,7 +262,7 @@ class LoopOverCoordinate(Node):
@staticmethod
def getLoopCounterSymbol(coordinateToLoopOver):
return TypedSymbol(LoopOverCoordinate.getLoopCounterName(coordinateToLoopOver), DataType('int'))
return TypedSymbol(LoopOverCoordinate.getLoopCounterName(coordinateToLoopOver), 'int')
@property
def loopCounterSymbol(self):
......@@ -479,8 +480,8 @@ class Indexed(Expr):
def __init__(self, args, base, parent=None):
super(Indexed, self).__init__(args, parent)
self.base = base
#Get dtype from label, and unpointer it
self.dtype = DataType(base.label.dtype.dtype)
# Get dtype from label, and unpointer it
self.dtype = createType(base.label.dtype.baseType)
def __repr__(self):
return '%s[%s]' % (self.args[0], self.args[1])
......
from .llvm import generateLLVM
from .cbackend import generateC, generateCUDA
try:
from .llvm import generateLLVM
except ImportError:
pass
from .cbackend import generateC
from .dot import dotprint
from sympy.utilities.codegen import CCodePrinter
from pystencils.astnodes import Node
from pystencils.types import createType
def generateC(astNode):
......@@ -7,17 +8,11 @@ def generateC(astNode):
Prints the abstract syntax tree as C function
"""
fieldTypes = set([f.dtype for f in astNode.fieldsAccessed])
useFloatConstants = "double" not in fieldTypes
printer = CBackend(cuda=False, constantsAsFloats=useFloatConstants)
useFloatConstants = createType("double") not in fieldTypes
printer = CBackend(constantsAsFloats=useFloatConstants)
return printer(astNode)
def generateCUDA(astNode):
fieldTypes = set([f.dtype for f in astNode.fieldsAccessed])
useFloatConstants = "double" not in fieldTypes
printer = CBackend(cuda=True, constantsAsFloats=useFloatConstants)
return printer(astNode)
# --------------------------------------- Backend Specific Nodes -------------------------------------------------------
......@@ -55,8 +50,7 @@ class PrintNode(CustomCppCode):
class CBackend(object):
def __init__(self, cuda=False, constantsAsFloats=False, sympyPrinter=None):
self.cuda = cuda
def __init__(self, constantsAsFloats=False, sympyPrinter=None):
if sympyPrinter is None:
self.sympyPrinter = CustomSympyPrinter(constantsAsFloats)
else:
......@@ -76,8 +70,7 @@ class CBackend(object):
def _print_KernelFunction(self, node):
functionArguments = ["%s %s" % (str(s.dtype), s.name) for s in node.parameters]
prefix = "__global__ void" if self.cuda else "void"
funcDeclaration = "%s %s(%s)" % (prefix, node.functionName, ", ".join(functionArguments))
funcDeclaration = "FUNC_PREFIX void %s(%s)" % (node.functionName, ", ".join(functionArguments))
body = self._print(node.body)
return funcDeclaration + "\n" + body
......
......@@ -6,7 +6,7 @@ from sympy import S
# S is numbers?
from pystencils.llvm.control_flow import Loop
from ..types import DataType
from ..types import createType, to_llvmlite_type
from ..astnodes import Indexed
......@@ -38,9 +38,9 @@ class LLVMPrinter(Printer):
self.tmp_var[name] = value
def _print_Number(self, n):
if n.dtype == DataType("int"):
if n.dtype == createType("int"):
return ir.Constant(self.integer, int(n))
elif n.dtype == DataType("double"):
elif n.dtype == createType("double"):
return ir.Constant(self.fp_type, float(n))
else:
raise NotImplementedError("Numbers can only have int and double", n)
......@@ -88,7 +88,7 @@ class LLVMPrinter(Printer):
def _print_Mul(self, expr):
nodes = [self._print(a) for a in expr.args]
e = nodes[0]
if expr.dtype == DataType('double'):
if expr.dtype == createType('double'):
mul = self.builder.fmul
else: # int TODO others?
mul = self.builder.mul
......@@ -99,7 +99,7 @@ class LLVMPrinter(Printer):
def _print_Add(self, expr):
nodes = [self._print(a) for a in expr.args]
e = nodes[0]
if expr.dtype == DataType('double'):
if expr.dtype == createType('double'):
add = self.builder.fadd
else: # int TODO others?
add = self.builder.add
......@@ -164,17 +164,22 @@ class LLVMPrinter(Printer):
from_dtype = conversion.args[0].dtype
# (From, to)
decision = {
(DataType("int"), DataType("double")): functools.partial(self.builder.sitofp, node, self.fp_type),
(DataType("double"), DataType("int")): functools.partial(self.builder.fptosi, node, self.integer),
(DataType("double *"), DataType("int")): functools.partial(self.builder.ptrtoint, node, self.integer),
(DataType("int"), DataType("double *")): functools.partial(self.builder.inttoptr, node, self.fp_pointer),
(DataType("double * __restrict__"), DataType("int")): functools.partial(self.builder.ptrtoint, node, self.integer),
(DataType("int"), DataType("double * __restrict__")): functools.partial(self.builder.inttoptr, node, self.fp_pointer),
(DataType("const double * __restrict__"), DataType("int")): functools.partial(self.builder.ptrtoint, node, self.integer),
(DataType("int"), DataType("const double * __restrict__")): functools.partial(self.builder.inttoptr, node, self.fp_pointer),
}
(createType("int"), createType("double")): functools.partial(self.builder.sitofp, node, self.fp_type),
(createType("double"), createType("int")): functools.partial(self.builder.fptosi, node, self.integer),
(createType("double *"), createType("int")): functools.partial(self.builder.ptrtoint, node, self.integer),
(createType("int"), createType("double *")): functools.partial(self.builder.inttoptr, node, self.fp_pointer),
(createType("double * restrict"), createType("int")): functools.partial(self.builder.ptrtoint, node, self.integer),
(createType("int"), createType("double * restrict")): functools.partial(self.builder.inttoptr, node, self.fp_pointer),
(createType("double * restrict const"), createType("int")): functools.partial(self.builder.ptrtoint, node, self.integer),
(createType("int"), createType("double * restrict const")): functools.partial(self.builder.inttoptr, node, self.fp_pointer),
}
# TODO float, const, restrict
# TODO bitcast, addrspacecast
# print([x for x in decision.keys()])
# print("Types:")
# print([(type(x), type(y)) for (x, y) in decision.keys()])
# print("Cast:")
# print((from_dtype, to_dtype))
return decision[(from_dtype, to_dtype)]()
def _print_Indexed(self, indexed):
......
This diff is collapsed.
import sympy as sp
from pystencils.transformations import resolveFieldAccesses, makeLoopOverDomain, typingFromSympyInspection, \
typeAllEquations, getOptimalLoopOrdering, parseBasePointerInfo, moveConstantsBeforeLoop, splitInnerLoop
from pystencils.types import TypedSymbol, DataType
from pystencils.types import TypedSymbol
from pystencils.field import Field
import pystencils.astnodes as ast
......@@ -37,7 +37,7 @@ def createKernel(listOfEquations, functionName="kernel", typeForSymbol=None, spl
if isinstance(term, Field.Access) or isinstance(term, TypedSymbol):
return term
elif isinstance(term, sp.Symbol):
return TypedSymbol(term.name, DataType(typeForSymbol[term.name]))
return TypedSymbol(term.name, typeForSymbol[term.name])
else:
raise ValueError("Term has to be field access or symbol")
......
import subprocess
import os
def getEnvironment(versionSpecifier, arch='x64'):
"""
Returns an environment dictionary, for activating the Visual Studio compiler
:param versionSpecifier: either a version number, year number, 'auto' or 'latest' for automatic detection of latest
installed version or 'setuptools' for setuptools-based detection
:param arch: x86 or x64
"""
if versionSpecifier == 'setuptools':
return getEnvironmentFromSetupTools(arch)
else:
if versionSpecifier in ('auto', 'latest'):
versionNr = findLatestMsvcVersionUsingEnvironmentVariables()
else:
versionNr = normalizeMsvcVersion(versionSpecifier)
vcVarsPath = getVcVarsPath(versionNr)
return getEnvironmentFromVcVarsFile(vcVarsPath, arch)
def findLatestMsvcVersionUsingEnvironmentVariables():
import re
regex = re.compile('VS(\d\d)\dCOMNTOOLS')
versions = []
for key, value in os.environ.items():
match = regex.match(key)
if match:
versions.append(int(match.group(1)))
if len(versions) == 0:
raise ValueError("Visual Studio not found.")
versions.sort()
return versions[-1]
def normalizeMsvcVersion(version):
"""
Takes version specifiers in the following form:
- year: 2012, 2013, 2015, either as int or string
- version numbers with or without dot i.e. 11.0 or 11
:return: integer version number
"""
if isinstance(version, str) and '.' in version:
version = version.split('.')[0]
version = int(version)
mapping = {
2015: 14,
2013: 12,
2012: 11
}
if version in mapping:
return mapping[version]
else:
return version
def getEnvironmentFromVcVarsFile(vcVarsFile, arch):
out = subprocess.check_output(
'cmd /u /c "{}" {} && set'.format(vcVarsFile, arch),
stderr=subprocess.STDOUT,
).decode('utf-16le', errors='replace')
env = {key.upper(): value for key, _, value in (line.partition('=') for line in out.splitlines()) if key and value}
return env
def getVcVarsPath(versionNr):
environmentVarName = 'VS%d0COMNTOOLS' % (versionNr,)
vcPath = os.environ[environmentVarName]
path = os.path.join(vcPath, '..', '..', 'VC', 'vcvarsall.bat')
return os.path.abspath(path)
def getEnvironmentFromSetupTools(arch):
from setuptools.msvc import msvc14_get_vc_env
msvcEnv = msvc14_get_vc_env(arch)
return {k.upper(): v for k, v in msvcEnv.items()}
......@@ -66,8 +66,8 @@ class EquationCollection(object):
newEquations = [fastSubs(eq, substitutionDict) for eq in self.mainEquations]
if addSubstitutionsAsSubexpressions:
newSubexpressions = [sp.Eq(b, a) for a, b in substitutionDict.items()] + newSubexpressions
return self.copy(newEquations, sortEquationsTopologically(newSubexpressions))
newSubexpressions = sortEquationsTopologically(newSubexpressions)
return self.copy(newEquations, newSubexpressions)
def addSimplificationHint(self, key, value):
"""
......
......@@ -3,7 +3,7 @@ import numpy as np
import sympy as sp
from sympy.core.cache import cacheit
from sympy.tensor import IndexedBase
from pystencils.types import TypedSymbol
from pystencils.types import TypedSymbol, createType
class Field(object):
......@@ -122,7 +122,7 @@ class Field(object):
def __init__(self, fieldName, dtype, layout, shape, strides):
"""Do not use directly. Use static create* methods"""
self._fieldName = fieldName
self._dtype = numpyDataTypeToC(dtype)
self._dtype = createType(dtype)
self._layout = normalizeLayout(layout)
self.shape = shape
self.strides = strides
......@@ -372,17 +372,6 @@ def computeStrides(shape, layout):
return tuple(strides)
def numpyDataTypeToC(dtype):
"""Mapping numpy data types to C data types"""
if dtype == np.float64:
return "double"
elif dtype == np.float32:
return "float"
elif dtype == np.int32:
return "int"
raise NotImplementedError("Cannot convert type " + str(dtype))
def offsetComponentToDirectionString(coordinateId, value):
"""
Translates numerical offset to string notation.
......
from pystencils.gpucuda.kernelcreation import createCUDAKernel
from pystencils.gpucuda.cudajit import makePythonFunction
from pystencils.backends.cbackend import generateCUDA
\ No newline at end of file
......@@ -3,7 +3,7 @@ import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
from pycuda.gpuarray import GPUArray
from pystencils.backends.cbackend import generateCUDA
from pystencils.backends.cbackend import generateC
from pystencils.transformations import symbolNameToVariableName
......@@ -58,7 +58,11 @@ def buildNumpyArgumentList(kernelFunctionNode, argumentDict):
def makePythonFunction(kernelFunctionNode, argumentDict={}):
mod = SourceModule(str(generateCUDA(kernelFunctionNode)), options=["-w"])
code = "#define FUNC_PREFIX __global__\n"
code += "#define RESTRICT __restrict__\n\n"
code += str(generateC(kernelFunctionNode))
mod = SourceModule(code, options=["-w"])
func = mod.get_function(kernelFunctionNode.functionName)
def wrapper(**kwargs):
......
import llvmlite.binding as llvm
import llvmlite.ir as ir
import llvmlite.binding as llvm
import logging.config
from ..types import toCtypes, createType
import ctypes as ct
......@@ -69,21 +71,19 @@ class Eval(object):
f.write(target_machine.emit_object(llvmmod))
fptr = {}
# TODO cpujit has its own string version
types = {str(ir.DoubleType()): ct.c_double,
str(ir.IntType(64)): ct.c_int64, # TODO int width?
str(ir.FloatType()): ct.c_float,
str(ir.VoidType()): None, # TODO thats a void pointer use None???
str(ir.DoubleType().as_pointer()): ct.POINTER(ct.c_double),
str(ir.IntType(64).as_pointer()): ct.POINTER(ct.c_int), # TODO int width?
str(ir.FloatType().as_pointer()): ct.POINTER(ct.c_float),
#ir.VoidType().as_pointer(): ct.c_void_p, # TODO there is no void pointer in llvmlite?
} # TODO Aggregate types
for function in module.functions:
if not function.is_declaration:
print(function.name)
print(type(function))
fptr[function.name] = ct.CFUNCTYPE(types[str(function.ftype.return_type)], *[types[str(arg)] for arg in function.ftype.args])(ee.get_function_address(function.name))
print(function.ftype.return_type)
print(type(function.ftype.return_type))
return_type = None
if function.ftype.return_type != ir.VoidType():
return_type = toCtypes(createType(str(function.ftype.return_type)))
args = [toCtypes(createType(str(arg))) for arg in function.ftype.args]
function_address = ee.get_function_address(function.name)
fptr[function.name] = ct.CFUNCTYPE(return_type, *args)(function_address)
# result = fptr(2, 3)
# print(result)
return fptr
......@@ -2,7 +2,7 @@ import sympy as sp
from pystencils.transformations import resolveFieldAccesses, makeLoopOverDomain, typingFromSympyInspection, \
typeAllEquations, getOptimalLoopOrdering, parseBasePointerInfo, moveConstantsBeforeLoop, splitInnerLoop, \
desympy_ast, insert_casts
from pystencils.types import TypedSymbol, DataType
from pystencils.types import TypedSymbol
from pystencils.field import Field
import pystencils.astnodes as ast
......@@ -36,7 +36,7 @@ def createKernel(listOfEquations, functionName="kernel", typeForSymbol=None, spl
if isinstance(term, Field.Access) or isinstance(term, TypedSymbol):
return term
elif isinstance(term, sp.Symbol):
return TypedSymbol(term.name, DataType(typeForSymbol[term.name]))
return TypedSymbol(term.name, typeForSymbol[term.name])
else:
raise ValueError("Term has to be field access or symbol")
......
......@@ -47,6 +47,10 @@ def normalizeSlice(slices, sizes):
return tuple(result)
def shiftSlice(slices, offset):
return [slice(k.start+offset, k.stop + offset, k.step) for k in slices]
def sliceFromDirection(directionName, dim, normalOffset=0, tangentialOffset=0):
"""
Create a slice from a direction named by compass scheme:
......@@ -98,3 +102,95 @@ def addGhostLayers(arr, indexDimensions=0, ghostLayers=1):
indexing += [slice(None, None, None)] * indexDimensions
result[indexing] = arr
return result
def getSliceBeforeGhostLayer(direction, ghostLayers=1, thickness=None, fullSlice=False):
"""
Returns slicing expression for region before ghost layer
:param direction: tuple specifying direction of slice
:param ghostLayers: number of ghost layers
:param thickness: thickness of the slice, defaults to number of ghost layers
:param fullSlice: if true also the ghost cells in directions orthogonal to direction are contained in the
returned slice. Example (d=W ): if fullSlice then also the ghost layer in N-S and T-B
are included, otherwise only inner cells are returned
"""
if not thickness:
thickness = ghostLayers
fullSliceInc = ghostLayers if not fullSlice else 0
slices = []
for dirComponent in direction:
if dirComponent == -1:
s = slice(ghostLayers, thickness + ghostLayers)
elif dirComponent == 0:
end = -fullSliceInc
s = slice(fullSliceInc, end if end != 0 else None)
elif dirComponent == 1:
start = -thickness - ghostLayers
end = -ghostLayers
s = slice(start if start != 0 else None, end if end != 0 else None)
else:
raise ValueError("Invalid direction: only -1, 0, 1 components are allowed")
slices.append(s)
return tuple(slices)
def getGhostRegionSlice(direction, ghostLayers=1, thickness=None, fullSlice=False):
"""
Returns slice of ghost region. For parameters see :func:`getSliceBeforeGhostLayer`
"""
if not thickness:
thickness = ghostLayers
assert thickness > 0
assert thickness <= ghostLayers
fullSliceInc = ghostLayers if not fullSlice else 0
slices = []
for dirComponent in direction:
if dirComponent == -1:
s = slice(ghostLayers - thickness, ghostLayers)
elif dirComponent == 0:
end = -fullSliceInc
s = slice(fullSliceInc, end if end != 0 else None)
elif dirComponent == 1:
start = -ghostLayers
end = - ghostLayers + thickness
s = slice(start if start != 0 else None, end if end != 0 else None)
else:
raise ValueError("Invalid direction: only -1, 0, 1 components are allowed")
slices.append(s)
return tuple(slices)
def getPeriodicBoundaryFunctor(stencil, ghostLayers=1, thickness=None):
"""
Returns a function that applies periodic boundary conditions
:param stencil: sequence of directions e.g. ( [0,1], [0,-1] ) for y periodicity
:param ghostLayers: how many ghost layers the array has
:param thickness: how many of the ghost layers to copy, None means 'all'
:return: function that takes a single array and applies the periodic copy operation
"""
srcDstSliceTuples = []
for d in stencil:
if sum([abs(e) for e in d]) == 0:
continue
invDir = (-e for e in d)
src = getSliceBeforeGhostLayer(invDir, ghostLayers, thickness=thickness, fullSlice=False)
dst = getGhostRegionSlice(d, ghostLayers, thickness=thickness, fullSlice=False)
print(d, ", ", src, "->", dst)
srcDstSliceTuples.append((src, dst))
def functor(field):
for srcSlice, dstSlice in srcDstSliceTuples:
field[dstSlice] = field[srcSlice]
return functor
if __name__ == '__main__':
import numpy as np
f = np.arange(0, 8*8).reshape(8, 8)
from lbmpy.stencils import getStencil
res = getPeriodicBoundaryFunctor(getStencil("D2Q9"), ghostLayers=2, thickness=1)
print(f)
res(f)
print(f)
......@@ -129,7 +129,7 @@ def replaceSecondOrderProducts(expr, searchSymbols, positive=None, replaceMixed=
else:
otherFactors *= t
if len(distinctVelTerms) == 2 and nrOfVelTerms == 2:
u, v = list(distinctVelTerms)
u, v = sorted(list(distinctVelTerms), key=lambda symbol: symbol.name)
if positive is None:
otherFactorsWithoutSymbols = otherFactors
for s in otherFactors.atoms(sp.Symbol):
......
from collections import defaultdict
from collections import defaultdict, OrderedDict
from operator import attrgetter
import sympy as sp
......@@ -6,7 +6,7 @@ from sympy.logic.boolalg import Boolean
from sympy.tensor import IndexedBase
from pystencils.field import Field, offsetComponentToDirectionString
from pystencils.types import TypedSymbol, DataType
from pystencils.types import TypedSymbol, createType, PointerType
from pystencils.slicing import normalizeSlice
import pystencils.astnodes as ast
......@@ -62,7 +62,7 @@ def makeLoopOverDomain(body, functionName, iterationSlice=None, ghostLayers=None
if len(shapeSet) != 1:
raise ValueError("Differently sized field accesses in loop body: " + str(shapeSet))
shape = list(shapeSet)[0]
shape = list(sorted(shapeSet, key=lambda e: str(e[0])))[0]
if iterationSlice is not None:
iterationSlice = normalizeSlice(iterationSlice, shape)
......@@ -109,7 +109,7 @@ def createIntermediateBasePointer(fieldAccess, coordinates, previousPtr):
Example:
>>> field = Field.createGeneric('myfield', spatialDimensions=2, indexDimensions=1)
>>> x, y = sp.symbols("x y")
>>> prevPointer = TypedSymbol("ptr", DataType("double"))
>>> prevPointer = TypedSymbol("ptr", "double")
>>> createIntermediateBasePointer(field[1,-2](5), {0: x}, prevPointer)
(ptr_E, x*fstride_myfield[0] + fstride_myfield[0])
>>> createIntermediateBasePointer(field[1,-2](5), {0: x, 1 : y }, prevPointer)
......@@ -140,7 +140,7 @@ def createIntermediateBasePointer(fieldAccess, coordinates, previousPtr):
if len(listToHash) > 0:
name += "%0.6X" % (abs(hash(tuple(listToHash))))
newPtr = TypedSymbol(previousPtr.name + name, DataType(previousPtr.dtype))
newPtr = TypedSymbol(previousPtr.name + name, previousPtr.dtype)
return newPtr, offset
......@@ -222,6 +222,9 @@ def resolveFieldAccesses(astNode, readOnlyFieldNames=set(), fieldToBasePointerIn
counters to index the field these symbols are used as coordinates
:return: transformed AST
"""
fieldToBasePointerInfo = OrderedDict(sorted(fieldToBasePointerInfo.items(), key=lambda pair: pair[0]))
fieldToFixedCoordinates = OrderedDict(sorted(fieldToFixedCoordinates.items(), key=lambda pair: pair[0]))
def visitSympyExpr(expr, enclosingBlock, sympyAssignment):
if isinstance(expr, Field.Access):
fieldAccess = expr
......@@ -231,12 +234,7 @@ def resolveFieldAccesses(astNode, readOnlyFieldNames=set(), fieldToBasePointerIn
else:
basePointerInfo = [list(range(field.indexDimensions + field.spatialDimensions))]
dtype = DataType(field.dtype)
dtype.alias = False
dtype.ptr = True
if field.name in readOnlyFieldNames:
dtype.const = True
dtype = PointerType(field.dtype, const=field.name in readOnlyFieldNames, restrict=True)
fieldPtr = TypedSymbol("%s%s" % (Field.DATA_PREFIX, symbolNameToVariableName(field.name)), dtype)
lastPointer = fieldPtr
......@@ -249,7 +247,7 @@ def resolveFieldAccesses(astNode, readOnlyFieldNames=set(), fieldToBasePointerIn
coordDict[e] = fieldToFixedCoordinates[field.name][e]
else:
ctrName = ast.LoopOverCoordinate.LOOP_COUNTER_NAME_PREFIX
coordDict[e] = TypedSymbol("%s_%d" % (ctrName, e), DataType('int'))
coordDict[e] = TypedSymbol("%s_%d" % (ctrName, e), 'int')
else:
coordDict[e] = fieldAccess.index[e-field.spatialDimensions]
return coordDict
......@@ -359,9 +357,8 @@ def splitInnerLoop(astNode, symbolGroups):
assert len(outerLoop) == 1, "Error in AST, multiple outermost loops."
outerLoop = outerLoop[0]
symbolsWithTemporaryArray = dict()
assignmentMap = {a.lhs: a for a in innerLoop.body.args}
symbolsWithTemporaryArray = OrderedDict()
assignmentMap = OrderedDict((a.lhs, a) for a in innerLoop.body.args)
assignmentGroups = []
for symbolGroup in symbolGroups:
......@@ -431,7 +428,7 @@ def typeAllEquations(eqs, typeForSymbol):
elif isinstance(term, TypedSymbol):
return term
elif isinstance(term, sp.Symbol):
return TypedSymbol(symbolNameToVariableName(term.name), DataType(typeForSymbol[term.name]))
return TypedSymbol(symbolNameToVariableName(term.name), typeForSymbol[term.name])
else:
newArgs = [processRhs(arg) for arg in term.args]
return term.func(*newArgs) if newArgs else term
......@@ -444,7 +441,7 @@ def typeAllEquations(eqs, typeForSymbol):
elif isinstance(term, TypedSymbol):
return term
elif isinstance(term, sp.Symbol):