ENH: add matrix-vector, vector-matrix multiplication (#1220)

- the vector-matrix multiplication is treated as a row vector

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -29,6 +29,81 @@ License
 #include <functional>
 #include <algorithm>
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Form, class Type>
+template<class ListType>
+Foam::tmp<Foam::Field<Type>> Foam::Matrix<Form, Type>::rightMultiplyImpl
+(
+    const ListType& colVec
+) const
+{
+    const Matrix<Form, Type>& mat = *this;
+
+    #ifdef FULLDEBUG
+    if (mat.n() != colVec.size())
+    {
+        FatalErrorInFunction
+            << "Attempt to multiply incompatible Matrix and Vector:" << nl
+            << "Matrix : (" << mat.m() << ", " << mat.n() << ')' << nl
+            << "Vector : " << colVec.size() << " rows" << nl
+            << "The number of Matrix columns must equal the Vector size" << nl
+            << abort(FatalError);
+    }
+    #endif
+
+    auto tresult = tmp<Field<Type>>::New(mat.m(), Zero);
+    auto& result = tresult.ref();
+
+    for (label i = 0; i < mat.m(); ++i)
+    {
+        for (label j = 0; j < mat.n(); ++j)
+        {
+            result[i] += mat(i, j)*colVec[j];
+        }
+    }
+
+    return tresult;
+}
+
+
+template<class Form, class Type>
+template<class ListType>
+Foam::tmp<Foam::Field<Type>> Foam::Matrix<Form, Type>::leftMultiplyImpl
+(
+    const ListType& rowVec
+) const
+{
+    const Matrix<Form, Type>& mat = *this;
+
+    #ifdef FULLDEBUG
+    if (rowVec.size() != mat.m())
+    {
+        FatalErrorInFunction
+            << "Attempt to multiply incompatible Matrix and Vector:" << nl
+            << "Matrix : (" << mat.m() << ", " << mat.n() << ')' << nl
+            << "Vector : " << rowVec.size() << " columns" << nl
+            << "The number of Matrix rows must equal the Vector size" << nl
+            << abort(FatalError);
+    }
+    #endif
+
+    auto tresult = tmp<Field<Type>>::New(mat.n(), Zero);
+    auto& result = tresult.ref();
+
+    for (label i = 0; i < mat.m(); ++i)
+    {
+        const Type& val = rowVec[i];
+        for (label j = 0; j < mat.n(); ++j)
+        {
+            result[j] += val*mat(i, j);
+        }
+    }
+
+    return tresult;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Form, class Type>
@@ -674,38 +749,6 @@ Foam::operator*
 }
 
 
-template<class Form, class Type>
-inline Foam::tmp<Foam::Field<Type>> Foam::operator*
-(
-    const Matrix<Form, Type>& mat,
-    const Field<Type>& x
-)
-{
-    if (mat.n() != x.size())
-    {
-        FatalErrorInFunction
-            << "Attempt to multiply incompatible matrix and field:" << nl
-            << "Matrix : (" << mat.m() << ", " << mat.n() << ')' << nl
-            << "Field  : " << x.size() << " rows" << nl
-            << "The number of matrix columns must equal the field size" << nl
-            << abort(FatalError);
-    }
-
-    auto tresult = tmp<Field<Type>>::New(mat.m(), Zero);
-    Field<Type>& result = tresult.ref();
-
-    for (label i=0; i < mat.m(); ++i)
-    {
-        for (label j=0; j < mat.n(); ++j)
-        {
-            result[i] += mat(i, j) * x[j];
-        }
-    }
-
-    return tresult;
-}
-
-
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "MatrixIO.C"

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -79,9 +79,20 @@ class Matrix
         //- Vector of values of type Type
         Type* __restrict__ v_;
 
+
+    // Private Member Functions
+
         //- Allocate storage for the contents
         inline void doAlloc();
 
+        //- Multiply matrix with a column vector on the right (A * x)
+        template<class ListType>
+        tmp<Field<Type>> rightMultiplyImpl(const ListType& colVec) const;
+
+        //- Multiply matrix with a row vector on the left (x * A)
+        template<class ListType>
+        tmp<Field<Type>> leftMultiplyImpl(const ListType& rowVec) const;
+
 
 public:
 
@@ -294,9 +305,37 @@ public:
         inline void shallowResize(const label m, const label n);
 
 
+    // Operations
+
         //- Return the transpose of the matrix
         Form T() const;
 
+        //- Multiply matrix with a column vector on the right (A * x)
+        inline tmp<Field<Type>> rightMultiply
+        (
+            const UList<Type>& colVec
+        ) const;
+
+        //- Multiply matrix with a column vector on the right (A * x)
+        template<class Addr>
+        inline tmp<Field<Type>> rightMultiply
+        (
+            const IndirectListBase<Type, Addr>& colVec
+        ) const;
+
+        //- Multiply matrix with a row vector on the left (x * A)
+        inline tmp<Field<Type>> leftMultiply
+        (
+            const UList<Type>& rowVec
+        ) const;
+
+        //- Multiply matrix with a row vector on the left (x * A)
+        template<class Addr>
+        inline tmp<Field<Type>> leftMultiply
+        (
+            const IndirectListBase<Type, Addr>& rowVec
+        ) const;
+
 
     // Member Operators
 
@@ -487,12 +526,37 @@ operator*
 );
 
 
-//- Matrix-vector multiplication
+//- Matrix-vector multiplication (A * x), where x is a column vector
 template<class Form, class Type>
-tmp<Field<Type>> operator*
+inline tmp<Field<Type>> operator*
 (
     const Matrix<Form, Type>& mat,
-    const Field<Type>& x
+    const UList<Type>& x
+);
+
+//- Matrix-vector multiplication (A * x), where x is a column vector
+template<class Form, class Type, class Addr>
+inline tmp<Field<Type>> operator*
+(
+    const Matrix<Form, Type>& mat,
+    const IndirectListBase<Type, Addr>& x
+);
+
+
+//- Vector-matrix multiplication (x * A), where x is a row vector
+template<class Form, class Type>
+inline tmp<Field<Type>> operator*
+(
+    const UList<Type>& x,
+    const Matrix<Form, Type>& mat
+);
+
+//- Vector-matrix multiplication (x * A), where x is a row vector
+template<class Form, class Type, class Addr>
+inline tmp<Field<Type>> operator*
+(
+    const IndirectListBase<Type, Addr>& x,
+    const Matrix<Form, Type>& mat
 );
 
 

src/OpenFOAM/matrices/Matrix/MatrixI.H

@@ -403,6 +403,48 @@ void Foam::Matrix<Form, Type>::shallowResize(const label m, const label n)
 }
 
 
+template<class Form, class Type>
+inline Foam::tmp<Foam::Field<Type>> Foam::Matrix<Form, Type>::rightMultiply
+(
+    const UList<Type>& colVec
+) const
+{
+    return this->rightMultiplyImpl(colVec);
+}
+
+
+template<class Form, class Type>
+template<class Addr>
+inline Foam::tmp<Foam::Field<Type>> Foam::Matrix<Form, Type>::rightMultiply
+(
+    const IndirectListBase<Type, Addr>& colVec
+) const
+{
+    return this->rightMultiplyImpl(colVec);
+}
+
+
+template<class Form, class Type>
+inline Foam::tmp<Foam::Field<Type>> Foam::Matrix<Form, Type>::leftMultiply
+(
+    const UList<Type>& rowVec
+) const
+{
+    return this->leftMultiplyImpl(rowVec);
+}
+
+
+template<class Form, class Type>
+template<class Addr>
+inline Foam::tmp<Foam::Field<Type>> Foam::Matrix<Form, Type>::leftMultiply
+(
+    const IndirectListBase<Type, Addr>& rowVec
+) const
+{
+    return this->leftMultiplyImpl(rowVec);
+}
+
+
 // * * * * * * * * * * * * * * * * * Iterators * * * * * * * * * * * * * * * //
 
 template<class Form, class Type>
@@ -497,4 +539,48 @@ inline Type* Foam::Matrix<Form, Type>::operator[](const label irow)
 }
 
 
+template<class Form, class Type>
+inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+(
+    const Matrix<Form, Type>& mat,
+    const UList<Type>& x
+)
+{
+    mat.rightMultiply(x);
+}
+
+
+template<class Form, class Type, class Addr>
+inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+(
+    const Matrix<Form, Type>& mat,
+    const IndirectListBase<Type, Addr>& x
+)
+{
+    mat.rightMultiply(x);
+}
+
+
+template<class Form, class Type>
+inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+(
+    const UList<Type>& x,
+    const Matrix<Form, Type>& mat
+)
+{
+    mat.leftMultiply(x);
+}
+
+
+template<class Form, class Type, class Addr>
+inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+(
+    const IndirectListBase<Type, Addr>& x,
+    const Matrix<Form, Type>& mat
+)
+{
+    mat.leftMultiply(x);
+}
+
+
 // ************************************************************************* //