Add transpose function

pyccel/ast/class_defs.py

@@ -10,7 +10,7 @@
 from .datatypes import (NativeBool, NativeInteger, NativeReal,
                         NativeComplex, NativeString)
 from .numpyext  import (Shape, NumpySum, NumpyAmin, NumpyAmax,
-                        NumpyImag, NumpyReal)
+                        NumpyImag, NumpyReal, NumpyTranspose)
 
 __all__ = ('BooleanClass',
         'IntegerClass',
@@ -135,6 +135,8 @@
         methods=[
             FunctionDef('shape',[],[],body=[],
                 decorators={'property':'property', 'numpy_wrapper':Shape}),
+            FunctionDef('T',[],[],body=[],
+                decorators={'property':'property', 'numpy_wrapper':NumpyTranspose}),
             FunctionDef('sum',[],[],body=[],
                 decorators={'numpy_wrapper':NumpySum}),
             FunctionDef('min',[],[],body=[],

pyccel/ast/numpyext.py

@@ -18,7 +18,7 @@
                              NativeReal, NativeComplex, NativeBool, str_dtype,
                              NativeNumeric)
 
-from .internals      import PyccelInternalFunction
+from .internals      import PyccelInternalFunction, Slice
 
 from .literals       import LiteralInteger, LiteralFloat, LiteralComplex, convert_to_literal
 from .literals       import LiteralTrue, LiteralFalse
@@ -950,7 +950,7 @@ class NumpyUfuncUnary(NumpyUfuncBase):
     def __init__(self, x):
         self._set_dtype_precision(x)
         self._set_shape_rank(x)
-        self._order      = x.order
+        self._set_order(x)
         super().__init__(x)
 
     def _set_shape_rank(self, x):
@@ -961,6 +961,9 @@ def _set_dtype_precision(self, x):
         self._dtype      = x.dtype if x.dtype is NativeComplex() else NativeReal()
         self._precision  = default_precision[str_dtype(self._dtype)]
 
+    def _set_order(self, x):
+        self._order      = x.order
+
 #------------------------------------------------------------------------------
 class NumpyUfuncBinary(NumpyUfuncBase):
     """Numpy's universal function with two arguments.
@@ -1173,6 +1176,46 @@ def _set_dtype_precision(self, x):
     def is_elemental(self):
         return False
 
+class NumpyTranspose(NumpyUfuncUnary):
+    """Represents a call to the transpose function in the Numpy library"""
+    __slots__ = ()
+    name = 'transpose'
+
+    def __new__(cls, x):
+        if x.rank<2:
+            return x
+        else:
+            return super().__new__(cls)
+
+    @property
+    def internal_var(self):
+        """ Return the variable being transposed
+        """
+        return self._args[0]
+
+    def __getitem__(self, *args):
+        x = self._args[0]
+        rank = x.rank
+        # Add empty slices to fully index the object
+        if len(args) < rank:
+            args = args + tuple([Slice(None, None)]*(rank-len(args)))
+        return NumpyTranspose(x.__getitem__(*reversed(args)))
+
+    def _set_dtype_precision(self, x):
+        self._dtype      = x.dtype
+        self._precision  = x.precision
+
+    def _set_shape_rank(self, x):
+        self._shape      = tuple(reversed(x.shape))
+        self._rank       = x.rank
+
+    def _set_order(self, x):
+        self._order = 'C' if x.order=='F' else 'F'
+
+    @property
+    def is_elemental(self):
+        return False
+
 #==============================================================================
 # TODO split numpy_functions into multiple dictionaries following
 # https://docs.scipy.org/doc/numpy-1.15.0/reference/routines.array-creation.html
@@ -1242,6 +1285,7 @@ def is_elemental(self):
     'arctanh'   : NumpyArctanh,
     # 'deg2rad'   : NumpyDeg2rad,
     # 'rad2deg'   : NumpyRad2deg,
+    'transpose' : NumpyTranspose,
 }
 
 numpy_linalg_functions = {

pyccel/ast/operators.py

@@ -183,8 +183,10 @@ def _set_order(self):
         Otherwise it defaults to 'C'
         """
         if self._rank is not None and self._rank > 1:
-            if all(a.order == self._args[0].order for a in self._args):
-                self._order = self._args[0].order
+            orders = [a.order for a in self._args if a.order is not None]
+            my_order = orders[0]
+            if all(o == my_order for o in orders):
+                self._order = my_order
             else:
                 self._order = 'C'
         else:

pyccel/ast/utilities.py

@@ -28,8 +28,9 @@
 from .literals      import LiteralString, LiteralInteger, Literal, Nil
 
 from .numpyext      import (NumpyEmpty, numpy_functions, numpy_linalg_functions,
-                            numpy_random_functions, numpy_constants, NumpyArray)
-from .operators     import PyccelAdd, PyccelMul, PyccelIs
+                            numpy_random_functions, numpy_constants, NumpyArray,
+                            NumpyTranspose)
+from .operators     import PyccelAdd, PyccelMul, PyccelIs, PyccelArithmeticOperator
 from .variable      import (Constant, Variable, ValuedVariable,
                             IndexedElement, InhomogeneousTupleVariable, VariableAddress,
                             HomogeneousTupleVariable )
@@ -233,7 +234,9 @@ def insert_index(expr, pos, index_var):
     >>> insert_index(expr, 0, i, language_has_vectors = True)
     c := a + b
     """
-    if isinstance(expr, (Variable, VariableAddress)):
+    if expr.rank==0:
+        return expr
+    elif isinstance(expr, (Variable, VariableAddress)):
         if expr.rank==0 or -pos>expr.rank:
             return expr
         if expr.shape[pos]==1:
@@ -244,6 +247,19 @@ def insert_index(expr, pos, index_var):
         indexes = [Slice(None,None)]*(expr.rank+pos) + [index_var]+[Slice(None,None)]*(-1-pos)
         return IndexedElement(expr, *indexes)
 
+    elif isinstance(expr, NumpyTranspose):
+        if expr.rank==0 or -pos>expr.rank:
+            return expr
+        if expr.shape[pos]==1:
+            # If there is no dimension in this axis, reduce the rank
+            index_var = LiteralInteger(0)
+
+        # Add index at the required position
+        if expr.rank<2:
+            return insert_index(expr.internal_var, expr.rank-1+pos, index_var)
+        else:
+            return NumpyTranspose(insert_index(expr.internal_var, expr.rank-1+pos, index_var))
+
     elif isinstance(expr, IndexedElement):
         base = expr.base
         indices = list(expr.indices)
@@ -268,6 +284,10 @@ def insert_index(expr, pos, index_var):
         indices[pos] = index_var
         return IndexedElement(base, *indices)
 
+    elif isinstance(expr, PyccelArithmeticOperator):
+        return type(expr)(insert_index(expr.args[0], pos, index_var),
+                          insert_index(expr.args[1], pos, index_var))
+
     else:
         raise NotImplementedError("Expansion not implemented for type : {}".format(type(expr)))
 
@@ -339,6 +359,7 @@ def collect_loops(block, indices, new_index_name, tmp_vars, language_has_vectors
             notable_nodes = line.get_attribute_nodes((Variable,
                                                        IndexedElement,
                                                        VariableAddress,
+                                                       NumpyTranspose,
                                                        FunctionCall,
                                                        PyccelInternalFunction,
                                                        PyccelIs))
@@ -357,22 +378,28 @@ def collect_loops(block, indices, new_index_name, tmp_vars, language_has_vectors
             variables      += [v for f in elemental_func_calls \
                                  for v in f.get_attribute_nodes((Variable, IndexedElement, VariableAddress),
                                                             excluded_nodes = (FunctionDef))]
+            transposed_vars = [v for v in notable_nodes if isinstance(v, NumpyTranspose)] \
+                                + [v for f in elemental_func_calls \
+                                     for v in f.get_attribute_nodes(NumpyTranspose,
+                                                            excluded_nodes = (FunctionDef))]
 
             is_checks = [n for n in notable_nodes if isinstance(n, PyccelIs)]
 
             variables = list(set(variables))
 
             # Check if the expression is already satisfactory
-            if compatible_operation(*variables, *is_checks, language_has_vectors = language_has_vectors):
+            if compatible_operation(*variables, *transposed_vars, *is_checks,
+                                    language_has_vectors = language_has_vectors):
                 result.append(line)
                 current_level = 0
                 continue
 
             # Find function calls in this line
-            funcs           = [f for f in notable_nodes if (isinstance(f, FunctionCall) \
+            funcs           = [f for f in notable_nodes+transposed_vars if (isinstance(f, FunctionCall) \
                                                             and not f.funcdef.is_elemental)]
-            internal_funcs  = [f for f in notable_nodes if (isinstance(f, PyccelInternalFunction) \
-                                                            and not f.is_elemental)]
+            internal_funcs  = [f for f in notable_nodes+transposed_vars if (isinstance(f, PyccelInternalFunction) \
+                                                            and not f.is_elemental) \
+                                                            and not isinstance(f, NumpyTranspose)]
 
             # Collect all variables for which values other than the value indexed in the loop are important
             # E.g. x = np.sum(a) has a dependence on a
@@ -390,8 +417,8 @@ def collect_loops(block, indices, new_index_name, tmp_vars, language_has_vectors
 
 
             if any(len(f.funcdef.results)!=1 for f in funcs):
-                errors.report("Loop unravelling cannot handle function calls \
-                        which return tuples or None",
+                errors.report("Loop unravelling cannot handle function calls "\
+                        "which return tuples or None",
                         symbol=line, severity='fatal')
 
             func_vars2 = [f.funcdef.results[0].clone(new_index_name('tmp')) for f in funcs]
@@ -400,9 +427,9 @@ def collect_loops(block, indices, new_index_name, tmp_vars, language_has_vectors
             if assigns:
                 # For now we do not handle memory allocation in loop unravelling
                 if any(v.rank > 0 for v in func_vars1) or any(v.rank > 0 for v in func_vars1):
-                    errors.report("Loop unravelling cannot handle extraction of function calls \
-                            which return arrays as this requires allocation. Please place the function \
-                            call on its own line",
+                    errors.report("Loop unravelling cannot handle extraction of function calls "\
+                            "which return arrays as this requires allocation. Please place the function "\
+                            "call on its own line",
                             symbol=line, severity='fatal')
                 line.substitute(internal_funcs, func_vars1, excluded_nodes=(FunctionCall))
                 line.substitute(funcs, func_vars2)
@@ -414,6 +441,7 @@ def collect_loops(block, indices, new_index_name, tmp_vars, language_has_vectors
             rank = line.lhs.rank
             shape = line.lhs.shape
             new_vars = variables
+            new_vars_t = transposed_vars
             # Loop over indexes, inserting until the expression can be evaluated
             # in the desired language
             new_level = 0
@@ -424,12 +452,15 @@ def collect_loops(block, indices, new_index_name, tmp_vars, language_has_vectors
                     indices.append(Variable('int',new_index_name('i')))
                 index_var = indices[rank+index]
                 new_vars = [insert_index(v, index, index_var) for v in new_vars]
-                if compatible_operation(*new_vars, language_has_vectors = language_has_vectors):
+                new_vars_t = [insert_index(v, index, index_var) for v in new_vars_t]
+                if compatible_operation(*new_vars, *new_vars_t, language_has_vectors = language_has_vectors):
                     break
 
             # Replace variable expressions with Indexed versions
             line.substitute(variables, new_vars, excluded_nodes = (FunctionCall, PyccelInternalFunction))
+            line.substitute(transposed_vars, new_vars_t, excluded_nodes = (FunctionCall))
             _ = [f.substitute(variables, new_vars, excluded_nodes = (FunctionDef)) for f in elemental_func_calls]
+            _ = [f.substitute(transposed_vars, new_vars_t, excluded_nodes = (FunctionDef)) for f in elemental_func_calls]
 
             # Recurse through result tree to save line with lines which need
             # the same set of for loops

pyccel/codegen/printing/ccode.py

@@ -1438,7 +1438,10 @@ def _print_AliasAssign(self, expr):
         # setting the pointer's is_view attribute to false so it can be ignored by the free_pointer function.
         if isinstance(lhs_var, Variable) and lhs_var.is_ndarray \
                 and isinstance(rhs_var, Variable) and rhs_var.is_ndarray:
-            return 'alias_assign(&{}, {});\n'.format(lhs, rhs)
+            if lhs_var.order == rhs_var.order:
+                return 'alias_assign(&{}, {});\n'.format(lhs, rhs)
+            else:
+                return 'transpose_alias_assign(&{}, {});\n'.format(lhs, rhs)
 
         return '{} = {};\n'.format(lhs, rhs)
 

pyccel/codegen/printing/fcode.py

@@ -2376,6 +2376,15 @@ def _print_NumpyUfuncBase(self, expr):
         code = '{0}({1})'.format(func_name, code_args)
         return self._get_statement(code)
 
+    def _print_NumpyTranspose(self, expr):
+        var = expr.internal_var
+        arg = self._print(var)
+        assign = expr.get_user_nodes(Assign)[0]
+        if assign.lhs.order != var.order:
+            return arg
+        else:
+            return 'transpose({0})'.format(arg)
+
     def _print_MathFunctionBase(self, expr):
         """ Convert a Python expression with a math function call to Fortran
         function call

pyccel/codegen/printing/pycode.py

@@ -13,7 +13,7 @@
 from pyccel.ast.datatypes  import default_precision
 from pyccel.ast.literals   import LiteralTrue, LiteralString
 from pyccel.ast.numpyext   import Shape as NumpyShape
-from pyccel.ast.variable   import DottedName, HomogeneousTupleVariable
+from pyccel.ast.variable   import DottedName, HomogeneousTupleVariable, Variable
 from pyccel.ast.utilities import builtin_import_registery as pyccel_builtin_import_registery
 
 from pyccel.codegen.printing.codeprinter import CodePrinter
@@ -38,6 +38,7 @@
                    'ones_like'  : 'ones',
                    'max'        : 'amax',
                    'min'        : 'amin',
+                   'T'          : 'transpose',
                    'full_like'  : 'full'},
         'numpy.random' : {'random' : 'rand'}
         }
@@ -475,14 +476,26 @@ def _print_Continue(self, expr):
         return 'continue\n'
 
     def _print_Assign(self, expr):
-        lhs = self._print(expr.lhs)
-        rhs = self._print(expr.rhs)
-        return'{0} = {1}\n'.format(lhs,rhs)
+        lhs = expr.lhs
+        rhs = expr.rhs
+
+        lhs_code = self._print(lhs)
+        rhs_code = self._print(rhs)
+        if isinstance(rhs, Variable) and lhs.rank>1 and rhs.order != lhs.order:
+            return'{0} = {1}.T\n'.format(lhs_code,rhs_code)
+        else:
+            return'{0} = {1}\n'.format(lhs_code,rhs_code)
 
     def _print_AliasAssign(self, expr):
-        lhs = self._print(expr.lhs)
-        rhs = self._print(expr.rhs)
-        return'{0} = {1}\n'.format(lhs,rhs)
+        lhs = expr.lhs
+        rhs = expr.rhs
+
+        lhs_code = self._print(lhs)
+        rhs_code = self._print(rhs)
+        if isinstance(rhs, Variable) and rhs.order!= lhs.order:
+            return'{0} = {1}.T\n'.format(lhs_code,rhs_code)
+        else:
+            return'{0} = {1}\n'.format(lhs_code,rhs_code)
 
     def _print_AugAssign(self, expr):
         lhs = self._print(expr.lhs)

pyccel/parser/semantic.py

@@ -81,6 +81,7 @@
 from pyccel.ast.numpyext import NumpyInt, NumpyInt8, NumpyInt16, NumpyInt32, NumpyInt64
 from pyccel.ast.numpyext import NumpyFloat, NumpyFloat32, NumpyFloat64
 from pyccel.ast.numpyext import NumpyComplex, NumpyComplex64, NumpyComplex128
+from pyccel.ast.numpyext import NumpyTranspose
 from pyccel.ast.numpyext import NumpyNewArray
 
 from pyccel.ast.omp import (OMP_For_Loop, OMP_Simd_Construct, OMP_Distribute_Construct,
@@ -771,6 +772,22 @@ def _infere_type(self, expr, **settings):
             d_var['cls_base'   ] = NumpyArrayClass
             return d_var
 
+        elif isinstance(expr, NumpyTranspose):
+
+            var = expr.internal_var
+
+            d_var['datatype'      ] = var.dtype
+            d_var['allocatable'   ] = var.allocatable
+            d_var['shape'         ] = tuple(reversed(var.shape))
+            d_var['rank'          ] = var.rank
+            d_var['cls_base'      ] = var.cls_base
+            d_var['is_pointer'    ] = var.is_pointer
+            d_var['is_target'     ] = var.is_target
+            d_var['order'         ] = 'C' if var.order=='F' else 'F'
+            d_var['precision'     ] = var.precision
+            d_var['is_stack_array'] = var.is_stack_array
+            return d_var
+
         elif isinstance(expr, PyccelAstNode):
 
             d_var['datatype'   ] = expr.dtype
@@ -1131,6 +1148,12 @@ def _ensure_target(self, rhs, d_lhs):
         """ Function using data about the new lhs to determine
         whether the lhs is a pointer and the rhs is a target
         """
+        if isinstance(rhs, NumpyTranspose) and rhs.internal_var.allocatable:
+            d_lhs['allocatable'] = False
+            d_lhs['is_pointer' ] = True
+            d_lhs['is_stack_array'] = False
+
+            rhs.internal_var.is_target = True
         if isinstance(rhs, Variable) and rhs.allocatable:
             d_lhs['allocatable'] = False
             d_lhs['is_pointer' ] = True
@@ -2217,6 +2240,10 @@ def _visit_Assign(self, expr, **settings):
                 d_var_i['shape'] = dvar['shape']
                 d_var_i['rank' ]  = dvar['rank']
 
+        elif isinstance(rhs, NumpyTranspose):
+            d_var  = self._infere_type(rhs, **settings)
+            rhs = rhs.internal_var
+
         else:
             d_var  = self._infere_type(rhs, **settings)
             d_list = d_var if isinstance(d_var, list) else [d_var]