diff --git a/orio/main/tuner/search/randomsimple/randomsimple.py b/orio/main/tuner/search/randomsimple/randomsimple.py
index 87c6fa12..471b07e9 100644
--- a/orio/main/tuner/search/randomsimple/randomsimple.py
+++ b/orio/main/tuner/search/randomsimple/randomsimple.py
@@ -14,7 +14,6 @@ def __init__(self, params):
         orio.main.tuner.search.search.Search.__init__(self, params)
 
         self.__readAlgoArgs()
-        #        self.have_z3 = False
 
         if self.time_limit <= 0 and self.total_runs <= 0:
             err((
@@ -24,9 +23,10 @@ def __init__(self, params):
     def searchBestCoord(self, startCoord=None):
 
         info('\n----- begin random search -----')
-        print("Total runs:", self.total_runs)
-        print("Time limit:", self.time_limit)
 
+        info( "Total runs: %d" % self.total_runs )
+        info( "Time limit: %d" % self.time_limit )
+        
         bestperfcost = self.MAXFLOAT
         bestcoord = None
         runs = 0
@@ -42,25 +42,28 @@ def searchBestCoord(self, startCoord=None):
             # get a random point
             coord = self.getRandomCoord()
 
-            # if not self.have_z3 and not self.checkValidity( coord ) or coord in visited:
-            if not self.checkValidity(coord) or coord in visited:
-                print("invalid point")
+            if coord == None:
+                debug( "No points left in the parameter space", obj=self, level=3 )
+                break
+            if not self.checkValidity( coord ) or coord in visited:
+                debug( "invalid point", obj=self, level=3 )
                 continue
             try:
-                print("coord:", coord, "run", runs)
-                perf_costs = self.getPerfCost(coord)
-                if bestperfcost > sum(perf_costs):
-                    info("Point %s gives a better perf: %s -- %s" % (coord, sum(perf_costs), bestperfcost))
-                    bestperfcost = sum(perf_costs)
+                debug( "coord: %s run %s" % (coord, runs ), obj=self, level=3 )
+                perf_costs = self.getPerfCost( coord )
+                if bestperfcost > sum( perf_costs ):
+                    info( "Point %s gives a better perf: %s -- %s" % (coord, sum( perf_costs ), bestperfcost ) )
+                    bestperfcost = sum( perf_costs )
                     bestcoord = coord
             except Exception as e:
                 info('FAILED: %s %s' % (e.__class__.__name__, e))
             runs += 1
-            # if not self.have_z3:
-            visited.append(coord)
-            # else:
-            #    point = self.coordToPerfParams( coord )
-            #    self.addPoint( point )
+
+            if not self.use_z3:
+                visited.append( coord )
+            else:
+                point = self.coordToPerfParams( coord )
+                self.z3solver.addPoint( point )
 
         search_time = time.time() - start_time
         return bestcoord, bestperfcost, search_time, runs
@@ -82,3 +85,4 @@ def __readAlgoArgs(self):
             else:
                 err('orio.main.tuner.search.randomsimple: unrecognized %s algorithm-specific argument: "%s"' %
                     (self.__class__.__name__, vname))
+
diff --git a/orio/main/tuner/search/search.py b/orio/main/tuner/search/search.py
index 75910dec..a4628b7e 100644
--- a/orio/main/tuner/search/search.py
+++ b/orio/main/tuner/search/search.py
@@ -5,7 +5,6 @@
 from orio.main.util.globals import *
 from functools import reduce
 
-
 class Search:
     '''The search engine used to explore the search space '''
 
@@ -75,12 +74,27 @@ def __init__(self, params):
         self.input_params = params.get('input_params')
         
         self.timing_code = ''
-        
+
         self.verbose = Globals().verbose
         self.perf_cost_records = {}
         self.transform_time={}
         self.best_coord_info="None"
 
+        # TODO pass it as an option
+        #        if 'use_z3' in params.keys():
+        try:
+            import z3_search
+            _use_z3 = True
+        except Exception as e:
+            _use_z3 = False
+
+        if _use_z3:
+            self.use_z3 = True
+            self.z3solver = z3_search.Z3search( self.total_dims, self.axis_names, self.axis_val_ranges, self.dim_uplimits, self.params['ptdriver'].tinfo.pparam_constraints, self )
+        else:
+            self.use_z3 = False
+            self.z3solver = None
+        
     #----------------------------------------------------------
 
     def searchBestCoord(self):
@@ -165,7 +179,17 @@ def coordToPerfParams(self, coord):
         perf_params = {}
         for i in range(0, self.total_dims):
             ipoint = coord[i]
-            perf_params[self.axis_names[i]] = self.axis_val_ranges[i][ipoint]
+            # If the point is not in the definition domain, take the nearest feasible value (z3)
+            if ipoint < len( self.axis_val_ranges[i] ):
+                perf_params[self.axis_names[i]] = self.axis_val_ranges[i][ipoint]
+            else:
+                if ipoint > self.axis_val_ranges[-1]:
+                    perf_params[self.axis_names[i]] = self.axis_val_ranges[i][-1]
+                else:
+                    perf_params[self.axis_names[i]] = self.axis_val_ranges[i][0]
+                    
+            # BN old version: perf_params[self.axis_names[i]] = self.axis_val_ranges[i][ipoint]
+
         return perf_params
 
     #----------------------------------------------------------
@@ -464,7 +488,13 @@ def getRandomCoord(self):
             iuplimit = self.dim_uplimits[i]
             ipoint = self.getRandomInt(0, iuplimit-1)
             random_coord.append(ipoint)
-        return random_coord
+        if self.use_z3:
+            point = self.z3solver.getNearestFeasible( random_coord )
+            if None == point:
+                return None
+            return self.z3solver.perfParamTabToCoord( point )
+        else:
+            return random_coord
                                                                      
 
     def getInitCoord(self):
@@ -475,9 +505,11 @@ def getInitCoord(self):
             #iuplimit = self.dim_uplimits[i]
             #ipoint = self.getRandomInt(0, iuplimit-1)
             random_coord.append(0)
-        return random_coord
-                                                                     
-
+        if self.use_z3:
+            point = self.z3solver.getNearestFeasible( random_coord )
+            return self.z3solver.perfParamTabToCoord( point )
+        else:
+            return random_coord
 
     #----------------------------------------------------------
 
diff --git a/orio/main/tuner/search/simplex/simplex.py b/orio/main/tuner/search/simplex/simplex.py
index e45734b3..17ec1a51 100644
--- a/orio/main/tuner/search/simplex/simplex.py
+++ b/orio/main/tuner/search/simplex/simplex.py
@@ -57,6 +57,8 @@ def __init__(self, params):
         self.x0 = [0] * self.total_dims
         self.sim_size = max(self.dim_uplimits)
 
+#        _have_z3 = False
+        
         # read all algorithm-specific arguments
         self.__readAlgoArgs()
 
@@ -430,7 +432,8 @@ def __readAlgoArgs(self):
                 err('man.tuner.search.simplex.simplex: unrecognized %s algorithm-specific argument: "%s"' %
                        (self.__class__.__name__, vname))
                 
-        
+    #-----------------------------------------------------
+
     #-----------------------------------------------------
 
     def __checkSearchSpace(self):
@@ -539,32 +542,68 @@ def __initRandomSimplex(self, simplex_records):
             
     #-----------------------------------------------------
 
+    # Get a centroid. If we are using z3, get a *feasible* centroid
+    # Get the nearest feasible point to the center
     def __getCentroid(self, coords):
-        '''Return a centroid coordinate'''
+        # Get a centroid
         total_coords = len(coords)
         centroid = coords[0]
         for c in coords[1:]:
             centroid = self.addCoords(centroid, c)
         centroid = self.mulCoords((1.0/total_coords), centroid)
-        return centroid
 
+        if self.use_z3:
+            # Get the nearest feasible point
+            point = self.z3solver.getNearestFeasible( centroid )
+            return self.z3solver.perfParamTabToCoord( point ) 
+        else:
+            return centroid
+
+    #-----------------------------------------------------
+    
     def __getReflection(self, coord, centroid):
         '''Return a reflection coordinate'''
         sub_coord = self.subCoords(centroid, coord)
-        return [self.addCoords(centroid, self.mulCoords(x, sub_coord)) for x in self.refl_coefs]
-    
+
+        point = map(lambda x: self.addCoords(centroid, self.mulCoords(x, sub_coord)),
+                    self.refl_coefs)
+        if self.use_z3:
+            # Get the nearest feasible point
+            return [ self.z3solver.perfParamTabToCoord( self.z3solver.getNearestFeasible( p ) ) for p in point ]
+       	else:
+            return list(point)
+
     def __getExpansion(self, coord, centroid):
         '''Return an expansion coordinate'''
         sub_coord = self.subCoords(coord, centroid)
-        return [self.addCoords(centroid, self.mulCoords(x, sub_coord)) for x in self.exp_coefs]
+        point = map(lambda x: self.addCoords(centroid, self.mulCoords(x, sub_coord)),
+                   self.exp_coefs)
+        if self.use_z3:
+            # Get the nearest feasible point
+            return [ self.z3solver.perfParamTabToCoord( self.z3solver.getNearestFeasible( p ) ) for p in point ]
+       	else:
+            return list(point)
     
     def __getContraction(self, coord, centroid):
         '''Return a contraction coordinate'''
         sub_coord = self.subCoords(coord, centroid)
-        return [self.addCoords(centroid, self.mulCoords(x, sub_coord)) for x in self.cont_coefs]
+
+        point = map(lambda x: self.addCoords(centroid, self.mulCoords(x, sub_coord)),
+                   self.cont_coefs)
+        if self.use_z3:
+            # Get the nearest feasible point
+            return [ self.z3solver.perfParamTabToCoord( self.z3solver.getNearestFeasible( p ) ) for p in point ]
+       	else:
+            return list(point)
 
     def __getShrinkage(self, coord, rest_coords):
         '''Return a shrinkage simplex'''
-        return [self.addCoords(coord, self.mulCoords(self.shri_coef,
-                                                                  self.subCoords(x, coord))) for x in rest_coords]
-    
+        point = map(lambda x: self.addCoords(coord, self.mulCoords(self.shri_coef,
+                                                                  self.subCoords(x, coord))),
+                   rest_coords)
+        if self.use_z3:
+            # Get the nearest feasible point
+            return [ self.z3solver.perfParamTabToCoord( self.z3solver.getNearestFeasible( p ) ) for p in point ]
+       	else:
+            return list(point)
+
diff --git a/orio/main/tuner/search/z3_search.py b/orio/main/tuner/search/z3_search.py
new file mode 100644
index 00000000..9bf8a741
--- /dev/null
+++ b/orio/main/tuner/search/z3_search.py
@@ -0,0 +1,305 @@
+from orio.main.util.globals import *
+
+try:
+    import z3
+    Globals.have_z3 = True
+except Exception as e:
+    Globals.have_z3 = False
+    raise ImportError( "Could not load z3. Will not be using it." )
+#    raise ImportError()  # Silent
+
+class Z3search:
+
+    def __init__( self,  total_dims, axis_names, axis_val_ranges, dim_uplimits, constraints, s ):
+        self.axis_names = axis_names
+        self.axis_val_ranges = axis_val_ranges
+        self.constraints = constraints
+        self.total_dims = total_dims
+        self.dim_uplimits = dim_uplimits
+        self.search = s
+
+        # Initialize the solver itself
+        self.variables = []
+        self.solver = z3.Solver()
+        self.optim = z3.Optimize()
+        Globals.z3types = {}
+        Globals.z3variables = {} 
+
+        # find what we have on the axis
+        self.__addVariableNames( self.axis_names, self.axis_val_ranges )
+        self.__addConstraints( self.constraints )
+
+    #-----------------------------------------------------
+        
+    # Defines the variables, their type and, for numeric ones, their definition domain
+    def __addVariableNames( self, names, ranges ):
+        for idx, name in enumerate( names ):
+            values = ranges[idx]
+            if ( [ True, False ] == values ) or ( [ False, True ] == values ): #len( values ) == 2 and ( False in values or True in values ):
+                Globals.z3variables[name] = z3.Bool( name )
+                self.__addDefinitionDomain( name, values )
+            else:
+                try:
+                    toto = int( values[0] )
+                    Globals.z3variables[name] = z3.Int( name )
+                    self.__addDefinitionDomain( name, values )
+                except ValueError:
+                    # if the parameter is non-numeric
+                    # TODO we can add a case to handle floats
+                    Globals.z3variables[name] = z3.Int( name )
+                    numvalues = self.__defineNonNumeric( name, values )
+                    self.__addDefinitionDomain( name, numvalues )
+            self.variables.append( name ) # FIXME variables is axis_names
+        return
+    
+    # Defines a variable's definition domain
+    def __addDefinitionDomain( self, var, dom ):
+        # get the variables
+        for k,v in Globals.z3variables.iteritems( ):
+            locals()[k] = v
+        if len( dom ) == 0:
+            return
+        definition = z3.Or( [ eval( "(" + var + " == " + str( v ) + ")" ) for v in dom ] )
+        # definition = z3.And( [ v >= 0, v < len( dom ) ] )
+        
+        self.solver.add( definition )
+        self.optim.add( definition )
+        return
+        
+    # Add the constraints
+    def __addConstraints( self, constraints ):
+        # get the variables
+        for k,v in Globals.z3variables.iteritems( ):
+            locals()[k] = v
+        for vname, rhs in constraints:
+            # convert to constraint to prefix syntax
+            toto = eval( self.__infixToPrefix( rhs ) )
+            self.solver.add( toto )
+            self.optim.add( toto )
+        return
+
+    # in case a parameter is not numeric (ie, alphanumeric), we cannot store it
+    # as a regular value with constraints in z3. Hence, put the values in a list and
+    # use their index in the list instead.
+    def __defineNonNumeric( self, name, values ):
+        Globals.z3types[name] = list( values )
+        return [ i for i in range( len( values ) ) ]
+    
+    # get the parameter corresponding to a numeric value
+    def __numericToNonNumeric( self, name, value, num ):
+        return Globals.z3types[name][num]
+    
+    # get the numeric value corresponding to a parameter
+    def __nonNumericToNumeric( self, name, value, param ):
+        return Globals.z3types[name].index( param )
+    
+    # Tests whether a z3 variable is numeric. Returns False for non-numeric ones.
+    def z3IsNumeric( self, dimensionNumber ):
+        return not self.axis_names[dimensionNumber] in Globals.z3types
+
+    def z3abs( self, x ):
+        return z3.If(x >= 0,x,-x)
+
+    #-----------------------------------------------------
+
+    # change constraint from infix notation (used in the input file) to prefix notation (used by z3)
+    def __infixToPrefix( self, expr ):
+        if "or" in expr or "and" in expr:
+            if "and" in expr:
+                operandsAND = expr.split( "and" )
+                return "z3.And(" + ", ".join( [ infixToPrefix( op ) if "or" in op or "and" in op else op  for op in operandsAND  ] ) + " ) "
+            elif "or" in expr:
+                operandsOR = expr.split( "or" )
+                return  "z3.Or(" + ", ".join( [ infixToPrefix( op ) if "or" in op or "and" in op else op for op in operandsOR  ] ) + " ) "
+            else:
+                return expr
+        else:
+            return expr
+        
+    #-----------------------------------------------------
+
+    # Add the constraints
+    def __addConstraints( self, constraints ):
+        # get the variables
+        for k,v in Globals.z3variables.iteritems( ):
+            locals()[k] = v
+        for vname, rhs in constraints:
+            # convert to constraint to prefix syntax
+            toto = eval( self.__infixToPrefix( rhs ) )
+            self.solver.add( toto )
+            self.optim.add( toto )
+        return
+
+    # Add a forbidden point in the solver (and the optimizer)
+    # Input: a dictionary of performance parameters.
+    def addPoint( self, point ):
+        # get the variables corresponding to the axis names
+        for name in self.axis_names:
+            locals()[name] = Globals.z3variables[name]
+        point2 = {}
+        # if there is a non-numeric value here, translate it
+        for idx, name in enumerate( point.keys() ):
+            # translate non-numeric values
+            if not self.z3IsNumeric( idx ): # non-numeric points -> coord to value
+                elem = Globals.z3types[name].index( point[name] )
+            else:
+                elem = point[name]
+            point2[name] = elem
+            
+        # remove this point
+        for s in [ self.solver, self.optim ]:
+            s.add( z3.Or( [ locals()[name] != point2[name] for name in self.axis_names ] ) )
+        return
+
+    #-----------------------------------------------------
+
+    # Returns a *feasible* random point
+    def getRandomCoord_z3( self ):
+        if self.solver.check() == z3.unsat:
+            return None
+        else:
+            model = self.solver.model()
+            return model
+
+    # Getting a random point from an SMT solver is not a trivial question.
+    # This one gives a somehow random point, but the distribution is not uniform.
+    # TODO: David's idea involving SAT.
+    # TODO: non-constrained variables?
+    def getRandomCoord_z3_distance( self ):
+        if self.solver.check() == z3.unsat:
+            return None
+        else:
+            # Pick a random point in the solution space
+            rpc = []
+            for i in range( self.total_dims ):
+                iuplimit = self.dim_uplimits[i]
+                ipoint = self.getRandomInt(0, iuplimit-1)
+                rpc.append(ipoint)
+
+            # Get the nearest feasible point
+            point = self.getNearestFeasible( rpc )
+            return self.perfParamTabToCoord( point )
+
+    # This is copied from search.py. Not exactly the cleanest way to do it, but simpler for the moment.
+    def getRandomInt(self, lbound, ubound):
+        '''To generate a random integer N such that lbound <= N <= ubound'''
+        from random import randint
+
+        if lbound > ubound:
+            err('orio.main.tuner.search.search internal error: the lower bound of genRandomInt must not be ' +
+                   'greater than the upper bound')
+        return randint(lbound, ubound)
+    
+    #-----------------------------------------------------
+
+    # For a given point, returns the nearest feasible point
+    def getNearestFeasible( self, coord ):
+        # get the variables corresponding to the axis names
+        for i in self.axis_names:
+            locals()[i] = Globals.z3variables[i]
+
+        # Convert into parameters
+        rpp = self.search.coordToPerfParams( coord )
+        for i, name in enumerate( self.axis_names ):
+            if not self.z3IsNumeric( i ):
+                rpp[name] = Globals.z3types[name].index( rpp[name] )
+                    
+        # create a new scope
+        self.optim.push()
+
+        # Get a possible point that minimizes the 1-norm distance to this random point
+        # forget the booleans
+        self.optim.minimize( z3.Sum( [ self.z3abs( locals()[name] - rpp[name] )for name in self.axis_names  if not z3.is_bool( locals()[name] ) ] ) )
+        if z3.unsat == self.optim.check():
+            return None
+        model = self.optim.model()
+
+        # restore the state (exit the scope)
+        self.optim.pop()
+        return self.z3ToPoint( model )
+
+    #-----------------------------------------------------
+        
+    # Defines the variables, their type and, for numeric ones, their definition domain
+    def __addVariableNames( self, names, ranges ):
+        for idx, name in enumerate( names ):
+            values = ranges[idx]
+            if ( [ True, False ] == values ) or ( [ False, True ] == values ): #len( values ) == 2 and ( False in values or True in values ):
+                Globals.z3variables[name] = z3.Bool( name )
+                self.__addDefinitionDomain( name, values )
+            else:
+                try:
+                    toto = int( values[0] )
+                    Globals.z3variables[name] = z3.Int( name )
+                    self.__addDefinitionDomain( name, values )
+                except ValueError:
+                    # if the parameter is non-numeric
+                    # TODO we can add a case to handle floats
+                    Globals.z3variables[name] = z3.Int( name )
+                    numvalues = self.__defineNonNumeric( name, values )
+                    self.__addDefinitionDomain( name, numvalues )
+            self.variables.append( name ) # FIXME variables is axis_names
+        return
+    
+    # z3-specific routine:
+    # Defines a variable's definition domain
+    def __addDefinitionDomain( self, var, dom ):
+        # get the variables
+        for k,v in Globals.z3variables.iteritems( ):
+            locals()[k] = v
+        if len( dom ) == 0:
+            return
+        definition = z3.Or( [ eval( "(" + var + " == " + str( v ) + ")" ) for v in dom ] )
+        # definition = z3.And( [ v >= 0, v < len( dom ) ] )
+        
+        self.solver.add( definition )
+        self.optim.add( definition )
+        return
+        
+    #-----------------------------------------------------
+
+    # Converts a model (returned by z3) into a point.
+    # Returns a list
+    def z3ToPoint( self, point ):
+        res = []
+        # get the variables corresponding to the axis names
+        for idx, i in enumerate( self.axis_names ):
+            locals()[i] = Globals.z3variables[i]
+            if i in Globals.z3types: # non-numeric points -> coord to value
+                value = Globals.z3types[i][point[locals()[i]].as_long()]
+            elif None == point[locals()[i]]:
+                info( "no value for %s, take %r (incomplete model)" % ( self.axis_names[idx], self.axis_val_ranges[idx][0] ) )
+                value = self.axis_val_ranges[idx][ 0]
+            elif z3.is_int_value( point[locals()[i]] ):
+                value = point[locals()[i]].as_long()
+            elif z3.is_true( point[locals()[i]] ) or z3.is_false( point[locals()[i]] ):
+                value = z3.is_true( point[locals()[i]] )
+            res.append( value )
+        return res
+
+    def perfParamToCoord( self, params ):
+        coord = [None]*self.total_dims
+        for i in range( self.total_dims ):
+            coord[i] = params[self.axis_names[i]]
+            if not self.z3IsNumeric( i ):
+                name = self.axis_names[i]
+                coord[i] = Globals.z3types[name].index( coord[i] )
+        return coord
+
+    def coordToTabOfPerfParams( self, coord ):
+        params = self.search.coordToPerfParams( coord )
+        perftab = []
+        for i, name in enumerate( self.axis_names ):
+            c = params[name]
+            if not self.z3IsNumeric( i ):
+                c = Globals.z3types[name].index( c )
+            perftab.append( c )
+        return perftab
+
+    def perfParamTabToCoord( self, params ):
+        coord = [None]*self.total_dims
+        for i in range( self.total_dims ):
+            coord[i] = self.axis_val_ranges[i].index( params[i] )
+        return coord
+
+