diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 4b3c44b8ebf1e8d091b766c38f71701435a0fcc9..589cbb781b5a3702ea625dcf9df59988c5fcee26 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/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"
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 700a31e523d303b805b2f1b30afdf08238e494d6..219630cd29d5285206ec438ffa376c049a070901 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/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 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 Field<Type>& x
+ 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
);
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index c264aed42ae5ff7a4d8df1f0d9e88e8c33b7369a..5c6d8e32221820d3c279747bca4ca7e067d22288 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/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);
+}
+
+
// ************************************************************************* //