[la] split VectorArrayInterface.dot into dot and pairwise_dot

src/pymor/algorithms/basisextension.py

@@ -174,7 +174,7 @@ def pod_basis_extension(basis, U, count=1, copy_basis=True, product=None, orthon
     new_basis = basis.copy() if copy_basis else basis
 
     if product is None:
-        U_proj_err = U - basis.lincomb(U.dot(basis, pairwise=False))
+        U_proj_err = U - basis.lincomb(U.dot(basis))
     else:
         U_proj_err = U - basis.lincomb(product.apply2(U, basis, pairwise=False))
 

src/pymor/algorithms/ei.py

@@ -130,7 +130,7 @@ def ei_greedy(U, error_norm=None, target_error=None, max_interpolation_dofs=None
         if projection == 'orthogonal':
             onb_collateral_basis.append(new_vec)
             gram_schmidt(onb_collateral_basis, offset=len(onb_collateral_basis) - 1, copy=False)
-            coeffs = ERR.dot(onb_collateral_basis, o_ind=len(onb_collateral_basis) - 1, pairwise=False)
+            coeffs = ERR.dot(onb_collateral_basis, o_ind=len(onb_collateral_basis) - 1)
             ERR.axpy(-coeffs[:, 0], onb_collateral_basis, x_ind=len(onb_collateral_basis) - 1)
 
     interpolation_matrix = collateral_basis.components(interpolation_dofs).T

src/pymor/la/genericsolvers.py

@@ -320,7 +320,7 @@ def lgmres(A, b, x0=None, tol=1e-5, maxiter=1000, M=None, callback=None,
             #     ++ orthogonalize
             hcur = []
             for v in vs:
-                alpha = v.dot(v_new, pairwise=False)[0, 0]
+                alpha = v.dot(v_new)[0, 0]
                 hcur.append(alpha)
                 v_new.axpy(-alpha, v)  # v_new -= alpha*v
             hcur.append(v_new.l2_norm()[0])

src/pymor/la/gram_schmidt.py

@@ -105,7 +105,7 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
                     if j in remove:
                         continue
                     if product is None:
-                        p = A.dot(A, ind=i, o_ind=j, pairwise=True)[0]
+                        p = A.pairwise_dot(A, ind=i, o_ind=j)[0]
                     else:
                         p = product.apply2(A, A, V_ind=i, U_ind=j, pairwise=True)[0]
                     A.axpy(-p, A, ind=i, x_ind=j)
@@ -131,7 +131,7 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
         if product:
             error_matrix = product.apply2(A, A, V_ind=range(offset, len(A)), pairwise=False)
         else:
-            error_matrix = A.dot(A, ind=range(offset, len(A)), pairwise=False)
+            error_matrix = A.dot(A, ind=range(offset, len(A)))
         error_matrix[:len(A) - offset, offset:len(A)] -= np.eye(len(A) - offset)
         if error_matrix.size > 0:
             err = np.max(np.abs(error_matrix))

src/pymor/la/interfaces.py

@@ -281,15 +281,13 @@ def axpy(self, alpha, x, ind=None, x_ind=None):
         pass
 
     @abstractmethod
-    def dot(self, other, pairwise, ind=None, o_ind=None):
+    def dot(self, other, ind=None, o_ind=None):
         """Returns the scalar products between |VectorArray| elements.
 
         Parameters
         ----------
         other
             A |VectorArray| containing the second factors.
-        pairwise
-            See return value documentation.
         ind
             Indices of the vectors whose scalar products are to be taken
             (see class documentation).
@@ -299,15 +297,34 @@ def dot(self, other, pairwise, ind=None, o_ind=None):
 
         Returns
         -------
-        If pairwise is `True`, returns a |NumPy array| `result` such
-        that ::
+        A |NumPy array| `result` such that ::
 
-            result[i] = ( self[ind][i], other[o_ind][i] ).
+            result[i, j] = ( self[ind][i], other[o_ind][j] ).
 
-        If pairwise is `False`, returns a |NumPy array| `result` such
-        that ::
+        """
+        pass
+
+    @abstractmethod
+    def pairwise_dot(self, other, ind=None, o_ind=None):
+        """Returns the scalar products between |VectorArray| elements.
+
+        Parameters
+        ----------
+        other
+            A |VectorArray| containing the second factors.
+        ind
+            Indices of the vectors whose scalar products are to be taken
+            (see class documentation).
+        o_ind
+            Indices of the vectors in `other` whose scalar products are to be
+            taken (see class documentation).
+
+        Returns
+        -------
+        A |NumPy array| `result` such that ::
+
+            result[i] = ( self[ind][i], other[o_ind][i] ).
 
-            result[i, j] = ( self[ind][i], other[o_ind][j] ).
         """
         pass
 
@@ -430,8 +447,8 @@ def amax(self, ind=None):
         pass
 
     def gramian(self, ind=None):
-        """Shorthand for `dot(self, pairwise=False, ind=ind, o_ind=ind)`."""
-        return self.dot(self, pairwise=False, ind=ind, o_ind=ind)
+        """Shorthand for `dot(self, ind=ind, o_ind=ind)`."""
+        return self.dot(self, ind=ind, o_ind=ind)
 
     def __add__(self, other):
         """The pairwise sum of two |VectorArrays|."""

src/pymor/la/listvectorarray.py

@@ -439,7 +439,7 @@ def axpy(self, alpha, x, ind=None, x_ind=None):
                 for xx, y in izip(X, Y):
                     y.axpy(alpha, xx)
 
-    def dot(self, other, pairwise, ind=None, o_ind=None):
+    def dot(self, other, ind=None, o_ind=None):
         assert self.check_ind(ind)
         assert other.check_ind(o_ind)
         assert self.space == other.space
@@ -464,17 +464,41 @@ def dot(self, other, pairwise, ind=None, o_ind=None):
             B = (other._list[i] for i in o_ind)
             len_B = len(o_ind)
 
-        if pairwise:
-            assert len_A == len_B
-            return np.array([a.dot(b) for a, b in izip(A, B)])
+        A = list(A)
+        B = list(B)
+        R = np.empty((len_A, len_B))
+        for i, a in enumerate(A):
+            for j, b in enumerate(B):
+                R[i, j] = a.dot(b)
+        return R
+
+    def pairwise_dot(self, other, ind=None, o_ind=None):
+        assert self.check_ind(ind)
+        assert other.check_ind(o_ind)
+        assert self.space == other.space
+
+        if ind is None:
+            A = self._list
+            len_A = len(A)
+        elif isinstance(ind, Number):
+            A = [self._list[ind]]
+            len_A = 1
         else:
-            A = list(A)
-            B = list(B)
-            R = np.empty((len_A, len_B))
-            for i, a in enumerate(A):
-                for j, b in enumerate(B):
-                    R[i, j] = a.dot(b)
-            return R
+            A = (self._list[i] for i in ind)
+            len_A = len(ind)
+
+        if o_ind is None:
+            B = other._list
+            len_B = len(B)
+        elif isinstance(o_ind, Number):
+            B = [other._list[o_ind]]
+            len_B = 1
+        else:
+            B = (other._list[i] for i in o_ind)
+            len_B = len(o_ind)
+
+        assert len_A == len_B
+        return np.array([a.dot(b) for a, b in izip(A, B)])
 
     def gramian(self, ind=None):
         assert self.check_ind(ind)

src/pymor/la/numpyvectorarray.py

@@ -254,7 +254,7 @@ def axpy(self, alpha, x, ind=None, x_ind=None):
             else:
                 self._array[ind] += (B * alpha)
 
-    def dot(self, other, pairwise, ind=None, o_ind=None):
+    def dot(self, other, ind=None, o_ind=None):
         assert self.check_ind(ind)
         assert other.check_ind(o_ind)
         assert self.dim == other.dim
@@ -270,11 +270,26 @@ def dot(self, other, pairwise, ind=None, o_ind=None):
         B = other._array[:other._len] if o_ind is None else \
             other._array[o_ind] if hasattr(o_ind, '__len__') else other._array[o_ind:o_ind + 1]
 
-        if pairwise:
-            assert self.len_ind(ind) == other.len_ind(o_ind)
-            return np.sum(A * B, axis=1)
-        else:
-            return A.dot(B.T)
+        return A.dot(B.T)
+
+    def pairwise_dot(self, other, ind=None, o_ind=None):
+        assert self.check_ind(ind)
+        assert other.check_ind(o_ind)
+        assert self.dim == other.dim
+        assert self.len_ind(ind) == other.len_ind(o_ind)
+
+        if NUMPY_INDEX_QUIRK:
+            if self._len == 0 and hasattr(ind, '__len__'):
+                ind = None
+            if other._len == 0 and hasattr(o_ind, '__len__'):
+                o_ind = None
+
+        A = self._array[:self._len] if ind is None else \
+            self._array[ind] if hasattr(ind, '__len__') else self._array[ind:ind + 1]
+        B = other._array[:other._len] if o_ind is None else \
+            other._array[o_ind] if hasattr(o_ind, '__len__') else other._array[o_ind:o_ind + 1]
+
+        return np.sum(A * B, axis=1)
 
     def lincomb(self, coefficients, ind=None):
         assert self.check_ind(ind)

src/pymor/la/pod.py

@@ -92,9 +92,9 @@ def pod(A, modes=None, product=None, tol=4e-8, symmetrize=False, orthonormalize=
         POD = gram_schmidt(POD, product=product, copy=False)
 
     if check:
-        if not product and not float_cmp_all(POD.dot(POD, pairwise=False), np.eye(len(POD)),
+        if not product and not float_cmp_all(POD.dot(POD), np.eye(len(POD)),
                                              atol=check_tol, rtol=0.):
-            err = np.max(np.abs(POD.dot(POD, pairwise=False) - np.eye(len(POD))))
+            err = np.max(np.abs(POD.dot(POD) - np.eye(len(POD))))
             raise AccuracyError('result not orthogonal (max err={})'.format(err))
         elif product and not float_cmp_all(product.apply2(POD, POD, pairwise=False), np.eye(len(POD)),
                                            atol=check_tol, rtol=0.):

src/pymor/operators/basic.py

@@ -37,7 +37,7 @@ def apply2(self, V, U, pairwise, U_ind=None, V_ind=None, mu=None, product=None):
         AU = self.apply(U, ind=U_ind, mu=mu)
         if product is not None:
             AU = product.apply(AU)
-        return V.dot(AU, ind=V_ind, pairwise=pairwise)
+        return V.pairwise_dot(AU, ind=V_ind) if pairwise else V.dot(AU, ind=V_ind)
 
     def jacobian(self, U, mu=None):
         if self.linear:

src/pymor/operators/constructions.py

@@ -346,7 +346,7 @@ def projected(self, source_basis, range_basis, product=None, name=None):
                 projected_value = NumpyVectorArray(product.apply2(range_basis, self._value, pairwise=False).T,
                                                    copy=False)
             else:
-                projected_value = NumpyVectorArray(range_basis.dot(self._value, pairwise=False).T, copy=False)
+                projected_value = NumpyVectorArray(range_basis.dot(self._value).T, copy=False)
         else:
             projected_value = self._value
         if source_basis is None:
@@ -447,7 +447,7 @@ def apply(self, U, ind=None, mu=None):
                 U = U.copy(ind)
             return self._array.lincomb(U.data)
         else:
-            return NumpyVectorArray(U.dot(self._array, ind=ind, pairwise=False), copy=False)
+            return NumpyVectorArray(U.dot(self._array, ind=ind), copy=False)
 
     def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None):
         assert U in self.range
@@ -457,7 +457,7 @@ def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product
             if range_product:
                 ATPrU = NumpyVectorArray(range_product.apply2(self._array, U, U_ind=ind, pairwise=False).T, copy=False)
             else:
-                ATPrU = NumpyVectorArray(self._array.dot(U, o_ind=ind, pairwise=False).T, copy=False)
+                ATPrU = NumpyVectorArray(self._array.dot(U, o_ind=ind).T, copy=False)
             if source_product:
                 return source_product.apply_inverse(ATPrU)
             else:

src/pymor/operators/ei.py

@@ -119,8 +119,7 @@ def projected(self, source_basis, range_basis, product=None, name=None):
 
             if range_basis is not None:
                 if product is None:
-                    projected_collateral_basis = NumpyVectorArray(self.collateral_basis.dot(range_basis,
-                                                                                            pairwise=False))
+                    projected_collateral_basis = NumpyVectorArray(self.collateral_basis.dot(range_basis))
                 else:
                     projected_collateral_basis = NumpyVectorArray(product.apply2(self.collateral_basis, range_basis,
                                                                                  pairwise=False))

src/pymor/playground/la/blockvectorarray.py

@@ -109,9 +109,19 @@ def axpy(self, alpha, x, ind=None, x_ind=None):
             for block, x_block in izip(self._blocks, x._blocks):
                 block.axpy(alpha, x_block, ind, x_ind)
 
-    def dot(self, other, pairwise, ind=None, o_ind=None):
+    def dot(self, other, ind=None, o_ind=None):
         assert other in self.space
-        dots = [block.dot(other_block, pairwise, ind=ind, o_ind=o_ind)
+        dots = [block.dot(other_block, ind=ind, o_ind=o_ind)
+                for block, other_block in izip(self._blocks, other._blocks)]
+        assert all([dot.shape == dots[0].shape for dot in dots])
+        ret = np.zeros(dots[0].shape)
+        for dot in dots:
+            ret += dot
+        return ret
+
+    def pairwise_dot(self, other, ind=None, o_ind=None):
+        assert other in self.space
+        dots = [block.pairwise_dot(other_block, ind=ind, o_ind=o_ind)
                 for block, other_block in izip(self._blocks, other._blocks)]
         assert all([dot.shape == dots[0].shape for dot in dots])
         ret = np.zeros(dots[0].shape)

src/pymor/playground/la/diskvectorarray.py

@@ -278,22 +278,27 @@ def axpy(self, alpha, x, ind=None, x_ind=None):
                     new_vec.axpy(alpha, x._load(xx))
                     self._store(y, new_vec)
 
-    def dot(self, other, pairwise, ind=None, o_ind=None):
+    def dot(self, other, ind=None, o_ind=None):
         assert self.check_ind(ind)
         assert other.check_ind(o_ind)
         assert self.space == other.space
         ind = list(xrange(self._len)) if ind is None else [ind] if isinstance(ind, Number) else ind
         o_ind = list(xrange(other._len)) if o_ind is None else [o_ind] if isinstance(o_ind, Number) else o_ind
 
-        if pairwise:
-            assert len(ind) == len(o_ind)
-            return np.array([self._load(i).dot(other._load(oi)) for i, oi in izip(ind, o_ind)])
-        else:
-            R = np.empty((len(ind), len(o_ind)))
-            for i, a in enumerate(ind):
-                for j, b in enumerate(o_ind):
-                    R[i, j] = self._load(a).dot(other._load(b))
-            return R
+        R = np.empty((len(ind), len(o_ind)))
+        for i, a in enumerate(ind):
+            for j, b in enumerate(o_ind):
+                R[i, j] = self._load(a).dot(other._load(b))
+        return R
+
+    def pairwise_dot(self, other, ind=None, o_ind=None):
+        assert self.check_ind(ind)
+        assert other.check_ind(o_ind)
+        assert self.space == other.space
+        ind = list(xrange(self._len)) if ind is None else [ind] if isinstance(ind, Number) else ind
+        o_ind = list(xrange(other._len)) if o_ind is None else [o_ind] if isinstance(o_ind, Number) else o_ind
+        assert len(ind) == len(o_ind)
+        return np.array([self._load(i).dot(other._load(oi)) for i, oi in izip(ind, o_ind)])
 
     def gramian(self, ind=None):
         assert self.check_ind(ind)

src/pymor/reductors/linear.py

@@ -131,9 +131,9 @@ def append_vector(U, R, RR):
                 append_vector(-op.apply(RB, [i]), R_O, RR_O)
 
     # compute Gram matrix of the residuals
-    R_RR = RR_R.dot(R_R, pairwise=False)
-    R_RO = np.hstack([RR_R.dot(R_O, pairwise=False) for R_O in R_Os])
-    R_OO = np.vstack([np.hstack([RR_O.dot(R_O, pairwise=False) for R_O in R_Os]) for RR_O in RR_Os])
+    R_RR = RR_R.dot(R_R)
+    R_RO = np.hstack([RR_R.dot(R_O) for R_O in R_Os])
+    R_OO = np.vstack([np.hstack([RR_O.dot(R_O) for R_O in R_Os]) for RR_O in RR_Os])
 
     estimator_matrix = np.empty((len(R_RR) + len(R_OO),) * 2)
     estimator_matrix[:len(R_RR), :len(R_RR)] = R_RR

src/pymor/reductors/residual.py

@@ -192,7 +192,7 @@ def apply(self, U, ind=None, mu=None):
         R = super(NonProjectedResiudalOperator, self).apply(U, ind=ind, mu=mu)
         if self.product:
             R_riesz = self.product.apply_inverse(R)
-            return R_riesz * (np.sqrt(R_riesz.dot(R, pairwise=True)) / R_riesz.l2_norm())[0]
+            return R_riesz * (np.sqrt(R_riesz.dot(R)) / R_riesz.l2_norm())[0]
         else:
             return R