diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C index c69bcc5ed2..de27a21058 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -29,6 +29,81 @@ License #include #include +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +template +Foam::tmp> Foam::Matrix::rightMultiplyImpl +( + const ListType& colVec +) const +{ + const Matrix& 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>::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 +template +Foam::tmp> Foam::Matrix::leftMultiplyImpl +( + const ListType& rowVec +) const +{ + const Matrix& 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>::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 @@ -674,38 +749,6 @@ Foam::operator* } -template -inline Foam::tmp> Foam::operator* -( - const Matrix& mat, - const Field& 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>::New(mat.m(), Zero); - Field& 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 9d5d709935..3ef51cbc84 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 + tmp> rightMultiplyImpl(const ListType& colVec) const; + + //- Multiply matrix with a row vector on the left (x * A) + template + tmp> 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> rightMultiply + ( + const UList& colVec + ) const; + + //- Multiply matrix with a column vector on the right (A * x) + template + inline tmp> rightMultiply + ( + const IndirectListBase& colVec + ) const; + + //- Multiply matrix with a row vector on the left (x * A) + inline tmp> leftMultiply + ( + const UList& rowVec + ) const; + + //- Multiply matrix with a row vector on the left (x * A) + template + inline tmp> leftMultiply + ( + const IndirectListBase& rowVec + ) const; + // Member Operators @@ -487,12 +526,37 @@ operator* ); -//- Matrix-vector multiplication +//- Matrix-vector multiplication (A * x), where x is a column vector template -tmp> operator* +inline tmp> operator* ( const Matrix& mat, - const Field& x + const UList& x +); + +//- Matrix-vector multiplication (A * x), where x is a column vector +template +inline tmp> operator* +( + const Matrix& mat, + const IndirectListBase& x +); + + +//- Vector-matrix multiplication (x * A), where x is a row vector +template +inline tmp> operator* +( + const UList& x, + const Matrix& mat +); + +//- Vector-matrix multiplication (x * A), where x is a row vector +template +inline tmp> operator* +( + const IndirectListBase& x, + const Matrix& mat ); diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H index 6b6519cb80..8abd40af82 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -403,6 +403,48 @@ void Foam::Matrix::shallowResize(const label m, const label n) } +template +inline Foam::tmp> Foam::Matrix::rightMultiply +( + const UList& colVec +) const +{ + return this->rightMultiplyImpl(colVec); +} + + +template +template +inline Foam::tmp> Foam::Matrix::rightMultiply +( + const IndirectListBase& colVec +) const +{ + return this->rightMultiplyImpl(colVec); +} + + +template +inline Foam::tmp> Foam::Matrix::leftMultiply +( + const UList& rowVec +) const +{ + return this->leftMultiplyImpl(rowVec); +} + + +template +template +inline Foam::tmp> Foam::Matrix::leftMultiply +( + const IndirectListBase& rowVec +) const +{ + return this->leftMultiplyImpl(rowVec); +} + + // * * * * * * * * * * * * * * * * * Iterators * * * * * * * * * * * * * * * // template @@ -497,4 +539,48 @@ inline Type* Foam::Matrix::operator[](const label irow) } +template +inline Foam::tmp> Foam::operator* +( + const Matrix& mat, + const UList& x +) +{ + mat.rightMultiply(x); +} + + +template +inline Foam::tmp> Foam::operator* +( + const Matrix& mat, + const IndirectListBase& x +) +{ + mat.rightMultiply(x); +} + + +template +inline Foam::tmp> Foam::operator* +( + const UList& x, + const Matrix& mat +) +{ + mat.leftMultiply(x); +} + + +template +inline Foam::tmp> Foam::operator* +( + const IndirectListBase& x, + const Matrix& mat +) +{ + mat.leftMultiply(x); +} + + // ************************************************************************* //