Type promotion

src/pymor/operators/block.py

@@ -139,9 +139,9 @@ def assemble_lincomb(self, operators, coefficients, solver_options=None, name=No
                     return None
             return self.__class__(blocks)
         else:
-            if coefficients[0] == 1:
+            c = coefficients[0]
+            if c == 1:
                 return self
-            c = coefficients[0].real if coefficients[0].imag == 0 else coefficients[0]
             for (i, j) in np.ndindex(self._blocks.shape):
                 blocks[i, j] = self._blocks[i, j] * c
             return self.__class__(blocks)
@@ -167,34 +167,24 @@ def __init__(self, blocks):
 
     def apply(self, U, mu=None):
         assert U in self.source
-
         V_blocks = [self._blocks[i, i].apply(U.block(i), mu=mu) for i in range(self.num_range_blocks)]
-
         return self.range.make_array(V_blocks)
 
     def apply_transpose(self, V, mu=None):
         assert V in self.range
-
         U_blocks = [self._blocks[i, i].apply_transpose(V.block(i), mu=mu) for i in range(self.num_source_blocks)]
-
-        U = self.source.make_array(U_blocks)
-
-        return U
+        return self.source.make_array(U_blocks)
 
     def apply_inverse(self, V, mu=None, least_squares=False):
         assert V in self.range
-
         U_blocks = [self._blocks[i, i].apply_inverse(V.block(i), mu=mu, least_squares=least_squares)
                     for i in range(self.num_source_blocks)]
-
         return self.source.make_array(U_blocks)
 
     def apply_inverse_transpose(self, U, mu=None, least_squares=False):
         assert U in self.source
-
         V_blocks = [self._blocks[i, i].apply_inverse_transpose(U.block(i), mu=mu, least_squares=least_squares)
                     for i in range(self.num_source_blocks)]
-
         return self.range.make_array(V_blocks)
 
     def assemble(self, mu=None):
@@ -224,9 +214,9 @@ def assemble_lincomb(self, operators, coefficients, solver_options=None, name=No
                     return None
             return self.__class__(blocks)
         else:
-            if coefficients[0] == 1:
+            c = coefficients[0]
+            if c == 1:
                 return self
-            c = coefficients[0].real if coefficients[0].imag == 0 else coefficients[0]
             for i in range(self.num_source_blocks):
                 blocks[i] = self._blocks[i, i] * c
             return self.__class__(blocks)

src/pymor/operators/constructions.py

@@ -5,6 +5,7 @@
 
 """Module containing some constructions to obtain new operators from old ones."""
 
+from functools import reduce
 from itertools import chain
 
 import numpy as np
@@ -86,18 +87,28 @@ def apply(self, U, mu=None):
 
     def apply2(self, V, U, mu=None):
         coeffs = self.evaluate_coefficients(mu)
-        R = self.operators[0].apply2(V, U, mu=mu)
-        R *= coeffs[0]
-        for op, c in zip(self.operators[1:], coeffs[1:]):
-            R += c * op.apply2(V, U, mu=mu)
+        matrices = [op.apply2(V, U, mu=mu) for op in self.operators]
+        coeffs_dtype = reduce(np.promote_types, (type(c) for c in coeffs))
+        matrices_dtype = reduce(np.promote_types, (m.dtype for m in matrices))
+        common_dtype = np.promote_types(coeffs_dtype, matrices_dtype)
+        R = coeffs[0] * matrices[0]
+        if R.dtype != common_dtype:
+            R = R.astype(common_dtype)
+        for m, c in zip(matrices[1:], coeffs[1:]):
+            R += c * m
         return R
 
     def pairwise_apply2(self, V, U, mu=None):
         coeffs = self.evaluate_coefficients(mu)
-        R = self.operators[0].pairwise_apply2(V, U, mu=mu)
-        R *= coeffs[0]
-        for op, c in zip(self.operators[1:], coeffs[1:]):
-            R += c * op.pairwise_apply2(V, U, mu=mu)
+        vectors = [op.pairwise_apply2(V, U, mu=mu) for op in self.operators]
+        coeffs_dtype = reduce(np.promote_types, (type(c) for c in coeffs))
+        vectors_dtype = reduce(np.promote_types, (v.dtype for v in vectors))
+        common_dtype = np.promote_types(coeffs_dtype, vectors_dtype)
+        R = coeffs[0] * vectors[0]
+        if R.dtype != common_dtype:
+            R = R.astype(common_dtype)
+        for v, c in zip(vectors[1:], coeffs[1:]):
+            R += c * v
         return R
 
     def apply_transpose(self, V, mu=None):
@@ -992,7 +1003,7 @@ def __init__(self, product, raise_negative, tol, name):
         self.build_parameter_type(product)
 
     def __call__(self, U, mu=None):
-        norm_squared = self.product.pairwise_apply2(U, U, mu=mu)
+        norm_squared = self.product.pairwise_apply2(U, U, mu=mu).real
         if self.tol > 0:
             norm_squared = np.where(np.logical_and(0 > norm_squared, norm_squared > - self.tol),
                                     0, norm_squared)

src/pymor/operators/numpy.py

@@ -337,16 +337,13 @@ def assemble_lincomb(self, operators, coefficients, solver_options=None, name=No
 
         common_mat_dtype = reduce(np.promote_types,
                                   (op._matrix.dtype for op in operators if hasattr(op, '_matrix')))
-        common_coef_dtype = reduce(np.promote_types, (type(c.real if c.imag == 0 else c) for c in coefficients))
+        common_coef_dtype = reduce(np.promote_types, (type(c) for c in coefficients))
         common_dtype = np.promote_types(common_mat_dtype, common_coef_dtype)
 
         if coefficients[0] == 1:
             matrix = operators[0]._matrix.astype(common_dtype)
         else:
-            if coefficients[0].imag == 0:
-                matrix = operators[0]._matrix * coefficients[0].real
-            else:
-                matrix = operators[0]._matrix * coefficients[0]
+            matrix = operators[0]._matrix * coefficients[0]
             if matrix.dtype != common_dtype:
                 matrix = matrix.astype(common_dtype)
 
@@ -371,11 +368,6 @@ def assemble_lincomb(self, operators, coefficients, solver_options=None, name=No
                     matrix -= op._matrix
                 except NotImplementedError:
                     matrix = matrix - op._matrix
-            elif c.imag == 0:
-                try:
-                    matrix += (op._matrix * c.real)
-                except NotImplementedError:
-                    matrix = matrix + (op._matrix * c.real)
             else:
                 try:
                     matrix += (op._matrix * c)

src/pymor/vectorarrays/block.py

@@ -3,6 +3,7 @@
 # Copyright 2013-2016 pyMOR developers and contributors. All rights reserved.
 # License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
 
+from functools import reduce
 from numbers import Number
 import numpy as np
 
@@ -30,6 +31,14 @@ def __init__(self, blocks, space):
     def data(self):
         return np.hstack([block.data for block in self._blocks])
 
+    @property
+    def real(self):
+        return BlockVectorArray([block.real for block in self._blocks], self.space)
+
+    @property
+    def imag(self):
+        return BlockVectorArray([block.imag for block in self._blocks], self.space)
+
     def block(self, ind):
         """
         Returns a copy of each block (no slicing).
@@ -84,7 +93,8 @@ def dot(self, other):
         assert other in self.space
         dots = [block.dot(other_block) for block, other_block in zip(self._blocks, other._blocks)]
         assert all([dot.shape == dots[0].shape for dot in dots])
-        ret = np.zeros(dots[0].shape)
+        common_dtype = reduce(np.promote_types, (dot.dtype for dot in dots))
+        ret = np.zeros(dots[0].shape, dtype=common_dtype)
         for dot in dots:
             ret += dot
         return ret
@@ -94,7 +104,8 @@ def pairwise_dot(self, other):
         dots = [block.pairwise_dot(other_block)
                 for block, other_block in zip(self._blocks, other._blocks)]
         assert all([dot.shape == dots[0].shape for dot in dots])
-        ret = np.zeros(dots[0].shape)
+        common_dtype = reduce(np.promote_types, (dot.dtype for dot in dots))
+        ret = np.zeros(dots[0].shape, dtype=common_dtype)
         for dot in dots:
             ret += dot
         return ret

src/pymor/vectorarrays/constructions.py

@@ -4,7 +4,7 @@
 
 
 def cat_arrays(vector_arrays):
-    """Return a new |VectorArray| which a concatenation of the arrays in `vector_arrays`."""
+    """Return a new |VectorArray| which is a concatenation of the arrays in `vector_arrays`."""
     vector_arrays = list(vector_arrays)
     total_length = sum(map(len, vector_arrays))
     cated_arrays = vector_arrays[0].empty(reserve=total_length)

src/pymor/vectorarrays/numpy.py

@@ -42,11 +42,11 @@ def data(self):
 
     @property
     def real(self):
-        return NumpyVectorArray(self._array[:self._len].real.copy(), self.space)
+        return NumpyVectorArray(self.data.real.copy(), self.space)
 
     @property
     def imag(self):
-        return NumpyVectorArray(self._array[:self._len].imag, self.space)
+        return NumpyVectorArray(self.data.imag.copy(), self.space)
 
     def __len__(self):
         return self._len
@@ -281,6 +281,10 @@ def __iadd__(self, other):
         assert self.dim == other.dim
         if self._refcount[0] > 1:
             self._deep_copy()
+        other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
+        common_dtype = np.promote_types(self._array.dtype, other_dtype)
+        if self._array.dtype != common_dtype:
+            self._array = self._array.astype(common_dtype)
         self._array[:self._len] += other.base._array[other.ind] if other.is_view else other._array[:other._len]
         return self
 
@@ -296,6 +300,10 @@ def __isub__(self, other):
         assert self.dim == other.dim
         if self._refcount[0] > 1:
             self._deep_copy()
+        other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
+        common_dtype = np.promote_types(self._array.dtype, other_dtype)
+        if self._array.dtype != common_dtype:
+            self._array = self._array.astype(common_dtype)
         self._array[:self._len] -= other.base._array[other.ind] if other.is_view else other._array[:other._len]
         return self
 
@@ -309,6 +317,10 @@ def __imul__(self, other):
             or isinstance(other, np.ndarray) and other.shape == (len(self),)
         if self._refcount[0] > 1:
             self._deep_copy()
+        other_dtype = other.dtype if isinstance(other, np.ndarray) else type(other)
+        common_dtype = np.promote_types(self._array.dtype, other_dtype)
+        if self._array.dtype != common_dtype:
+            self._array = self._array.astype(common_dtype)
         self._array[:self._len] *= other
         return self
 
@@ -463,7 +475,6 @@ def l2_norm(self):
 
     def l2_norm2(self):
         return self.base.l2_norm2(_ind=self.ind)
-        pass
 
     def sup_norm(self):
         return self.base.sup_norm(_ind=self.ind)
@@ -488,6 +499,10 @@ def __iadd__(self, other):
         assert self.base.check_ind_unique(self.ind)
         if self.base._refcount[0] > 1:
             self._deep_copy()
+        other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
+        common_dtype = np.promote_types(self.base._array.dtype, other_dtype)
+        if self.base._array.dtype != common_dtype:
+            self.base._array = self.base._array.astype(common_dtype)
         self.base.array[self.ind] += other.base._array[other.ind] if other.is_view else other._array[:other._len]
         return self
 
@@ -504,6 +519,10 @@ def __isub__(self, other):
         assert self.base.check_ind_unique(self.ind)
         if self.base._refcount[0] > 1:
             self._deep_copy()
+        other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
+        common_dtype = np.promote_types(self.base._array.dtype, other_dtype)
+        if self.base._array.dtype != common_dtype:
+            self.base._array = self.base._array.astype(common_dtype)
         self.base._array[self.ind] -= other.base._array[other.ind] if other.is_view else other._array[:other._len]
         return self
 
@@ -518,6 +537,10 @@ def __imul__(self, other):
         assert self.base.check_ind_unique(self.ind)
         if self.base._refcount[0] > 1:
             self._deep_copy()
+        other_dtype = other.dtype if isinstance(other, np.ndarray) else type(other)
+        common_dtype = np.promote_types(self.base._array.dtype, other_dtype)
+        if self.base._array.dtype != common_dtype:
+            self.base._array = self.base._array.astype(common_dtype)
         self.base._array[self.ind] *= other
         return self
 

src/pymortests/complex_values.py

@@ -24,7 +24,7 @@ def test_complex():
     # assemble_lincomb
     assert not np.iscomplexobj(Aop.assemble_lincomb((Iop, Bop), (1, 1))._matrix)
     assert not np.iscomplexobj(Aop.assemble_lincomb((Aop, Bop), (1, 1))._matrix)
-    assert not np.iscomplexobj(Aop.assemble_lincomb((Aop, Bop), (1 + 0j, 1 + 0j))._matrix)
+    assert np.iscomplexobj(Aop.assemble_lincomb((Aop, Bop), (1 + 0j, 1 + 0j))._matrix)
     assert np.iscomplexobj(Aop.assemble_lincomb((Aop, Bop), (1j, 1))._matrix)
     assert np.iscomplexobj(Aop.assemble_lincomb((Bop, Aop), (1, 1j))._matrix)