From 061eb53fb5a1c6420289629324a1185b81877772 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Wed, 22 May 2019 12:18:31 +0100 Subject: [PATCH] ENH: improvements, modernization of matrix containers (#1220) - add iterators, begin/end, empty() methods for STL behaviour. Use standard algorithms where possible * std::fill, std::copy * std::min_element, std::max_element - access methods consistent with other OpenFOAM containers: * data(), cdata(), uniform() - Use ListPolicy to impose output line breaks - Can recover matrix storage for re-use elsewhere. For example, to populate values with 2D i-j addressing and later release it as flat linear storage. - construct/assign moveable - added minMax() function for Matrix - additional inplace +=, -=, *=, /= operations - add named methods at() and rowData() to Matrix. Allows a better distinction between linear and row-based addressing - low-level matrix solve on List/UList instead of Field --- applications/test/Matrix/Test-Matrix.C | 114 ++-- .../matrices/DiagonalMatrix/DiagonalMatrix.C | 76 ++- .../matrices/DiagonalMatrix/DiagonalMatrix.H | 33 +- src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C | 40 +- src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H | 19 +- .../matrices/LUscalarMatrix/LUscalarMatrix.H | 4 +- .../LUscalarMatrix/LUscalarMatrixTemplates.C | 37 +- src/OpenFOAM/matrices/Matrix/Matrix.C | 585 ++++++++++-------- src/OpenFOAM/matrices/Matrix/Matrix.H | 341 ++++++---- src/OpenFOAM/matrices/Matrix/MatrixI.H | 215 ++++++- src/OpenFOAM/matrices/Matrix/MatrixIO.C | 176 +++--- src/OpenFOAM/matrices/QRMatrix/QRMatrix.C | 31 +- src/OpenFOAM/matrices/QRMatrix/QRMatrix.H | 25 +- src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H | 2 +- .../RectangularMatrix/RectangularMatrix.H | 15 +- .../RectangularMatrix/RectangularMatrixI.H | 15 +- .../matrices/SquareMatrix/SquareMatrix.H | 38 +- .../matrices/SquareMatrix/SquareMatrixI.H | 49 +- .../SymmetricSquareMatrix.H | 12 +- .../SymmetricSquareMatrixI.H | 24 +- 20 files changed, 1121 insertions(+), 730 deletions(-) diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C index 7c9cb61195..a402c73a30 100644 --- a/applications/test/Matrix/Test-Matrix.C +++ b/applications/test/Matrix/Test-Matrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -40,44 +40,82 @@ using namespace Foam; int main(int argc, char *argv[]) { - SquareMatrix hmm(3); - - hmm(0, 0) = -3.0; - hmm(0, 1) = 10.0; - hmm(0, 2) = -4.0; - hmm(1, 0) = 2.0; - hmm(1, 1) = 3.0; - hmm(1, 2) = 10.0; - hmm(2, 0) = 2.0; - hmm(2, 1) = 6.0; - hmm(2, 2) = 1.0; - - //Info<< hmm << endl << hmm - 2.0*(-hmm) << endl; - Info<< max(hmm) << endl; - Info<< min(hmm) << endl; - - SquareMatrix hmm2(3, I); - - hmm = hmm2; - - Info<< hmm << endl; - - //SquareMatrix hmm3(Sin); - - //Info<< hmm3 << endl; - - SquareMatrix hmm4; - - hmm4 = hmm2; - - Info<< hmm4 << endl; - - SquareMatrix hmm5; - - hmm4 = hmm5; - Info<< hmm5 << endl; - + // SquareMatrix { + Info<< nl << "Test SquareMatrix" << nl; + + SquareMatrix square1(3); + + square1(0, 0) = -3.0; + square1(0, 1) = 10.0; + square1(0, 2) = -4.0; + square1(1, 0) = 2.0; + square1(1, 1) = 3.0; + square1(1, 2) = 10.0; + square1(2, 0) = 2.0; + square1(2, 1) = 6.0; + square1(2, 2) = 1.0; + + Info<< "matrix: " << square1 + << " begin: " << long(square1.cdata()) << nl; + + //Info<< square1 - 2.0*(-square1) << nl; + Info<< "min:" << min(square1) << " max:" << max(square1) << nl; + Info<< "min/max: " << minMax(square1) << nl; + + // Steal matrix contents + + List stole(square1.release()); + + Info<< "matrix: " << square1 + << " begin: " << long(square1.cdata()) << nl; + + Info<< "List: " << stole << nl; + + Info<< "min/max: " << minMax(square1) << nl; + + square1 = 100; + + Info<< "matrix: " << square1 + << " begin: " << long(square1.cdata()) << nl; + + + SquareMatrix square2(3, I); + + square1 = square2; + + Info<< nl << "After copy assign from Identity:" << nl << square1 << nl; + + square1 += 1.25*square2; + + Info<< nl << "After +=" << nl << square1 << nl; + + square1 -= square2.T(); + + Info<< nl << "After -=" << nl << square1 << nl; + + square1 *= 10; + + Info<< nl << "After *=" << nl << square1 << nl; + + square1 /= 8; + + Info<< nl << "After /=" << nl << square1 << nl; + + SquareMatrix square4; + square4 = square2; + + Info<< nl << square4 << endl; + + SquareMatrix square5; + square4 = square5; + Info<< nl << square5 << endl; + } + + // RectangularMatrix + { + Info<< nl << "Test RectangularMatrix" << nl; + RectangularMatrix rm1(5, 6, 3.1); rm1(0, 1) = 4.5; RectangularMatrix rm1b(rm1.block(2, 2, 0, 0)); diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C index 6e4dad4b56..f05aede750 100644 --- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -37,45 +37,54 @@ inline Foam::DiagonalMatrix::DiagonalMatrix() template -template -Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) +Foam::DiagonalMatrix::DiagonalMatrix(const label n) : - List(min(a.m(), a.n())) + List(n) +{} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label n, const zero) +: + List(n, Zero) +{} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label n, const Type& val) +: + List(n, val) +{} + + +template +template +Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& mat) +: + List(min(mat.m(), mat.n())) { - forAll(*this, i) + label i = 0; + + for (Type& val : *this) { - this->operator[](i) = a(i, i); + val = mat(i, i); + ++i; } } -template -Foam::DiagonalMatrix::DiagonalMatrix(const label size) -: - List(size) -{} - - -template -Foam::DiagonalMatrix::DiagonalMatrix(const label size, const Type& val) -: - List(size, val) -{} - - template Foam::DiagonalMatrix& Foam::DiagonalMatrix::invert() { - forAll(*this, i) + for (Type& val : *this) { - Type x = this->operator[](i); - if (mag(x) < VSMALL) + if (mag(val) < VSMALL) { - this->operator[](i) = Type(0); + val = Zero; } else { - this->operator[](i) = Type(1)/x; + val = Type(1)/val; } } @@ -84,21 +93,24 @@ Foam::DiagonalMatrix& Foam::DiagonalMatrix::invert() template -Foam::DiagonalMatrix Foam::inv(const DiagonalMatrix& A) +Foam::DiagonalMatrix Foam::inv(const DiagonalMatrix& mat) { - DiagonalMatrix Ainv = A; + DiagonalMatrix Ainv(mat.size()); - forAll(A, i) + Type* iter = Ainv.begin(); + + for (const Type& val : mat) { - Type x = A[i]; - if (mag(x) < VSMALL) + if (mag(val) < VSMALL) { - Ainv[i] = Type(0); + *iter = Zero; } else { - Ainv[i] = Type(1)/x; + *iter = Type(1)/val; } + + ++iter; } return Ainv; diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H index 2521693d1f..a7ba834b9a 100644 --- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -27,8 +27,7 @@ Class Foam::DiagonalMatrix Description - DiagonalMatrix is a 2D diagonal matrix of objects - of type Type, size nxn + A 2D diagonal matrix of objects of type \, size (N x N) SourceFiles DiagonalMatrix.C @@ -45,12 +44,11 @@ SourceFiles namespace Foam { -// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * // - +// Forward Declarations template class Matrix; /*---------------------------------------------------------------------------*\ - Class DiagonalMatrix Declaration + Class DiagonalMatrix Declaration \*---------------------------------------------------------------------------*/ template @@ -62,24 +60,29 @@ public: // Constructors - //- Null constructor. + //- Construct null. DiagonalMatrix(); - //- Construct from diagonal component of a Matrix - template - DiagonalMatrix(const Matrix&); - //- Construct empty from size - DiagonalMatrix(const label size); + explicit DiagonalMatrix(const label n); + + //- Construct from size + //- initializing all elements to the given value + DiagonalMatrix(const label n, const zero); //- Construct from size and a value - DiagonalMatrix(const label, const Type&); + DiagonalMatrix(const label n, const Type& val); + + //- Construct from the diagonal of a Matrix + template + DiagonalMatrix(const Matrix& mat); - // Member functions + // Member Functions //- Invert the diagonal matrix and return itself DiagonalMatrix& invert(); + }; @@ -87,7 +90,7 @@ public: //- Return the diagonal Matrix inverse template -DiagonalMatrix inv(const DiagonalMatrix&); +DiagonalMatrix inv(const DiagonalMatrix& mat); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C index c18bbe5c51..7798d2ea9d 100644 --- a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C +++ b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -35,27 +35,27 @@ Foam::LLTMatrix::LLTMatrix() template -Foam::LLTMatrix::LLTMatrix(const SquareMatrix& M) +Foam::LLTMatrix::LLTMatrix(const SquareMatrix& mat) { - decompose(M); + decompose(mat); } // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // template -void Foam::LLTMatrix::decompose(const SquareMatrix& M) +void Foam::LLTMatrix::decompose(const SquareMatrix& mat) { SquareMatrix& LLT = *this; // Initialize the LLT decomposition matrix to M - LLT = M; + LLT = mat; const label m = LLT.m(); - for (label i=0; i i) { @@ -65,7 +65,7 @@ void Foam::LLTMatrix::decompose(const SquareMatrix& M) Type sum = LLT(i, j); - for (label k=0; k::decompose(const SquareMatrix& M) template void Foam::LLTMatrix::solve ( - Field& x, - const Field& source + List& x, + const UList& source ) const { - // If x and source are different initialize x = source + // If x and source are different, copy initialize x = source if (&x != &source) { x = source; @@ -106,11 +106,11 @@ void Foam::LLTMatrix::solve const SquareMatrix& LLT = *this; const label m = LLT.m(); - for (label i=0; i::solve x[i] = sum/LLT(i, i); } - for (int i=m - 1; i >= 0; i--) + for (label i=m - 1; i >= 0; --i) { Type sum = x[i]; - for (label j=i + 1; j Foam::tmp> Foam::LLTMatrix::solve ( - const Field& source + const UList& source ) const { - tmp> tx(new Field(this->m())); - Field& x = tx.ref(); + auto tresult(tmp>::New(source.size())); - solve(x, source); + solve(tresult.ref(), source); - return tx; + return tresult; } diff --git a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H index 80d42c23b2..a9f853e41c 100644 --- a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H +++ b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -58,7 +58,6 @@ class LLTMatrix : public SquareMatrix { - public: // Constructors @@ -66,23 +65,23 @@ public: //- Construct null LLTMatrix(); - //- Construct from a square matrix and perform the decomposition - LLTMatrix(const SquareMatrix& M); + //- Construct and perform the decomposition on input square matrix + LLTMatrix(const SquareMatrix& mat); // Member Functions - //- Perform the Cholesky decomposition of the matrix - void decompose(const SquareMatrix& M); + //- Copy matrix and perform Cholesky decomposition + void decompose(const SquareMatrix& mat); //- Solve the linear system with the given source - // and returning the solution in the Field argument x. + //- and returning the solution in the Field argument x. // This function may be called with the same field for x and source. - void solve(Field& x, const Field& source) const; + void solve(List& x, const UList& source) const; //- Solve the linear system with the given source - // returning the solution - tmp> solve(const Field& source) const; + //- returning the solution + tmp> solve(const UList& source) const; }; diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H index 006ec01f6b..f03b4025f2 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H @@ -122,12 +122,12 @@ public: // and returning the solution in the Field argument x. // This function may be called with the same field for x and source. template - void solve(Field& x, const Field& source) const; + void solve(List& x, const UList& source) const; //- Solve the linear system with the given source // returning the solution template - tmp> solve(const Field& source) const; + tmp> solve(const UList& source) const; //- Set M to the inverse of this square matrix void inv(scalarSquareMatrix& M) const; diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C index f9329df0c7..584a86593a 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2017 OpenFOAM Foundation @@ -26,15 +26,15 @@ License \*---------------------------------------------------------------------------*/ #include "LUscalarMatrix.H" -#include "SubField.H" +#include "SubList.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template void Foam::LUscalarMatrix::solve ( - Field& x, - const Field& source + List& x, + const UList& source ) const { // If x and source are different initialize x = source @@ -45,15 +45,13 @@ void Foam::LUscalarMatrix::solve if (Pstream::parRun()) { - Field X(m()); + List X; // scratch space (on master) if (Pstream::master(comm_)) { - typename Field::subField - ( - X, - x.size() - ) = x; + X.resize(m()); + + SubList(X, x.size()) = x; for ( @@ -82,7 +80,7 @@ void Foam::LUscalarMatrix::solve ( Pstream::commsTypes::scheduled, Pstream::masterNo(), - reinterpret_cast(x.begin()), + reinterpret_cast(x.cdata()), x.byteSize(), Pstream::msgType(), comm_ @@ -93,11 +91,7 @@ void Foam::LUscalarMatrix::solve { LUBacksubstitute(*this, pivotIndices_, X); - x = typename Field::subField - ( - X, - x.size() - ); + x = SubList(X, x.size()); for ( @@ -114,7 +108,7 @@ void Foam::LUscalarMatrix::solve ( &(X[procOffsets_[slave]]) ), - (procOffsets_[slave + 1]-procOffsets_[slave])*sizeof(Type), + (procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type), Pstream::msgType(), comm_ ); @@ -126,7 +120,7 @@ void Foam::LUscalarMatrix::solve ( Pstream::commsTypes::scheduled, Pstream::masterNo(), - reinterpret_cast(x.begin()), + reinterpret_cast(x.data()), x.byteSize(), Pstream::msgType(), comm_ @@ -143,13 +137,12 @@ void Foam::LUscalarMatrix::solve template Foam::tmp> Foam::LUscalarMatrix::solve ( - const Field& source + const UList& source ) const { - tmp> tx(new Field(m())); - Field& x = tx.ref(); + auto tx(tmp>::New(m())); - solve(x, source); + solve(tx.ref(), source); return tx; } diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C index d9dbc580bd..c69bcc5ed2 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2017 OpenFOAM Foundation @@ -26,18 +26,8 @@ License \*---------------------------------------------------------------------------*/ #include "Matrix.H" - -// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // - -template -void Foam::Matrix::allocate() -{ - if (mRows_ && nCols_) - { - v_ = new Type[size()]; - } -} - +#include +#include // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -48,14 +38,9 @@ Foam::Matrix::Matrix(const label m, const label n) nCols_(n), v_(nullptr) { - if (mRows_ < 0 || nCols_ < 0) - { - FatalErrorInFunction - << "Incorrect m, n " << mRows_ << ", " << nCols_ - << abort(FatalError); - } + checkSize(); - allocate(); + doAlloc(); } @@ -66,90 +51,71 @@ Foam::Matrix::Matrix(const label m, const label n, const zero) nCols_(n), v_(nullptr) { - if (mRows_ < 0 || nCols_ < 0) - { - FatalErrorInFunction - << "Incorrect m, n " << mRows_ << ", " << nCols_ - << abort(FatalError); - } + checkSize(); - allocate(); + doAlloc(); - if (v_) - { - const label mn = size(); - for (label i=0; i -Foam::Matrix::Matrix(const label m, const label n, const Type& s) +Foam::Matrix::Matrix(const label m, const label n, const Type& val) : mRows_(m), nCols_(n), v_(nullptr) { - if (mRows_ < 0 || nCols_ < 0) - { - FatalErrorInFunction - << "Incorrect m, n " << mRows_ << ", " << nCols_ - << abort(FatalError); - } + checkSize(); - allocate(); + doAlloc(); - if (v_) + std::fill(begin(), end(), val); +} + + +template +Foam::Matrix::Matrix(const Matrix& mat) +: + mRows_(mat.mRows_), + nCols_(mat.nCols_), + v_(nullptr) +{ + if (mat.cdata()) { - const label mn = size(); - for (label i=0; i -Foam::Matrix::Matrix(const Matrix& M) +Foam::Matrix::Matrix(Matrix&& mat) : - mRows_(M.mRows_), - nCols_(M.nCols_), - v_(nullptr) + mRows_(mat.mRows_), + nCols_(mat.nCols_), + v_(mat.v_) { - if (M.v_) - { - allocate(); - - const label mn = size(); - for (label i=0; i template -Foam::Matrix::Matrix(const Matrix& M) +Foam::Matrix::Matrix(const Matrix& mat) : - mRows_(M.m()), - nCols_(M.n()), + mRows_(mat.m()), + nCols_(mat.n()), v_(nullptr) { - if (M.v()) + if (mat.cdata()) { - allocate(); + doAlloc(); - const label mn = size(); - for (label i=0; i::Matrix mRows_(Mb.m()), nCols_(Mb.n()) { - allocate(); + doAlloc(); - for (label i=0; i::Matrix mRows_(Mb.m()), nCols_(Mb.n()) { - allocate(); + doAlloc(); - for (label i=0; i::clear() template -void Foam::Matrix::transfer(Matrix& M) +Foam::List Foam::Matrix::release() { + List list; + + const label len = size(); + + if (v_ && len) + { + UList storage(v_, len); + list.swap(storage); + + v_ = nullptr; + } clear(); - mRows_ = M.mRows_; - M.mRows_ = 0; - - nCols_ = M.nCols_; - M.nCols_ = 0; - - v_ = M.v_; - M.v_ = nullptr; + return list; } template -void Foam::Matrix::setSize(const label m, const label n) +void Foam::Matrix::swap(Matrix& mat) { - mType newMatrix(m, n, Zero); + Foam::Swap(mRows_, mat.mRows_); + Foam::Swap(nCols_, mat.nCols_); + Foam::Swap(v_, mat.v_); +} - label minM = min(m, mRows_); - label minN = min(n, nCols_); - for (label i=0; i +void Foam::Matrix::transfer(Matrix& mat) +{ + clear(); + + mRows_ = mat.mRows_; + nCols_ = mat.nCols_; + v_ = mat.v_; + + mat.mRows_ = 0; + mat.nCols_ = 0; + mat.v_ = nullptr; +} + + +template +void Foam::Matrix::resize(const label m, const label n) +{ + if (m == mRows_ && n == nCols_) { - for (label j=0; j newMatrix(m, n, Zero); + + const label mrow = min(m, mRows_); + const label ncol = min(n, nCols_); + + for (label i=0; i < mrow; ++i) + { + for (label j=0; j < ncol; ++j) { newMatrix(i, j) = (*this)(i, j); } @@ -265,14 +264,13 @@ void Foam::Matrix::setSize(const label m, const label n) template Form Foam::Matrix::T() const { - const Matrix& A = *this; Form At(n(), m()); - for (label i=0; i::T() const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template -void Foam::Matrix::operator=(const Matrix& M) +void Foam::Matrix::operator=(const Matrix& mat) { - if (this == &M) + if (this == &mat) { FatalErrorInFunction << "Attempted assignment to self" << abort(FatalError); } - if (mRows_ != M.mRows_ || nCols_ != M.nCols_) + if (mRows_ != mat.mRows_ || nCols_ != mat.nCols_) { clear(); - mRows_ = M.mRows_; - nCols_ = M.nCols_; - allocate(); + mRows_ = mat.mRows_; + nCols_ = mat.nCols_; + doAlloc(); } if (v_) { - const label mn = size(); - for (label i=0; i +void Foam::Matrix::operator=(Matrix&& mat) +{ + if (this == &mat) + { + FatalErrorInFunction + << "Attempted assignment to self" + << abort(FatalError); + } + + this->transfer(mat); +} + + template template void Foam::Matrix::operator= @@ -318,11 +326,11 @@ void Foam::Matrix::operator= const ConstMatrixBlock& Mb ) { - for (label i=0; i::operator= const MatrixBlock& Mb ) { - for (label i=0; i -void Foam::Matrix::operator=(const Type& s) +void Foam::Matrix::operator=(const Type& val) { - if (v_) - { - const label mn = size(); - for (label i=0; i void Foam::Matrix::operator=(const zero) { - if (v_) + std::fill(begin(), end(), Zero); +} + + +template +void Foam::Matrix::operator+=(const Matrix& other) +{ + if (this == &other) { - const label mn = size(); - for (label i=0; idata(); + const Type* in = other.cdata(); + + const label len = this->size(); + + for (label idx=0; idx < len; ++idx) + { + out[idx] += in[idx]; + } +} + + +template +void Foam::Matrix::operator-=(const Matrix& other) +{ + if (this == &other) + { + FatalErrorInFunction + << "Attempted subtraction from self" + << abort(FatalError); + } + + if (m() != other.m() || n() != other.n()) + { + FatalErrorInFunction + << "Attempt to subtract matrices with different sizes: (" + << m() << ", " << n() << ") != (" + << other.m() << ", " << other.n() << ')' << nl + << abort(FatalError); + } + + Type* out = this->data(); + const Type* in = other.cdata(); + + const label len = this->size(); + + for (label idx=0; idx < len; ++idx) + { + out[idx] -= in[idx]; + } +} + + +template +void Foam::Matrix::operator*=(const scalar s) +{ + for (Type& val : *this) + { + val *= s; + } +} + + +template +void Foam::Matrix::operator/=(const scalar s) +{ + for (Type& val : *this) + { + val /= s; } } @@ -376,89 +452,61 @@ void Foam::Matrix::operator=(const zero) // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // template -const Type& Foam::max(const Matrix& M) +const Type& Foam::max(const Matrix& mat) { - const label mn = M.size(); - - if (mn) - { - label curMaxI = 0; - const Type* Mv = M.v(); - - for (label i=1; i Mv[curMaxI]) - { - curMaxI = i; - } - } - - return Mv[curMaxI]; - } - else + if (mat.empty()) { FatalErrorInFunction - << "Matrix is empty" - << abort(FatalError); - - // Return in error to keep compiler happy - return M(0, 0); + << "Matrix is empty" << abort(FatalError); } + + return *(std::max_element(mat.cbegin(), mat.cend())); } template -const Type& Foam::min(const Matrix& M) +const Type& Foam::min(const Matrix& mat) { - const label mn = M.size(); - - if (mn) - { - label curMinI = 0; - const Type* Mv = M.v(); - - for (label i=1; i +Foam::MinMax Foam::minMax(const Matrix& mat) +{ + MinMax result; + + for (const Type& val : mat) + { + result += val; + } + + return result; } // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // template -Form Foam::operator-(const Matrix& M) +Form Foam::operator-(const Matrix& mat) { - Form nM(M.m(), M.n()); + Form result(mat.m(), mat.n()); - if (M.m() && M.n()) - { - Type* nMv = nM.v(); - const Type* Mv = M.v(); + std::transform + ( + mat.cbegin(), + mat.cend(), + result.begin(), + std::negate() + ); - const label mn = M.size(); - for (label i=0; i& A, const Matrix& B) if (A.m() != B.m()) { FatalErrorInFunction - << "Attempt to add matrices with different numbers of rows: " - << A.m() << ", " << B.m() - << abort(FatalError); - } - - if (A.n() != B.n()) - { - FatalErrorInFunction - << "Attempt to add matrices with different numbers of columns: " - << A.n() << ", " << B.n() + << "Attempt to add matrices with different sizes: (" + << A.m() << ", " << A.n() << ") != (" + << B.m() << ", " << B.n() << ')' << nl << abort(FatalError); } Form AB(A.m(), A.n()); - Type* ABv = AB.v(); - const Type* Av = A.v(); - const Type* Bv = B.v(); + Type* ABv = AB.data(); + const Type* Av = A.cdata(); + const Type* Bv = B.cdata(); - const label mn = A.size(); - for (label i=0; i& A, const Matrix& B) if (A.m() != B.m()) { FatalErrorInFunction - << "Attempt to add matrices with different numbers of rows: " - << A.m() << ", " << B.m() - << abort(FatalError); - } - - if (A.n() != B.n()) - { - FatalErrorInFunction - << "Attempt to add matrices with different numbers of columns: " - << A.n() << ", " << B.n() + << "Attempt to subtract matrices with different sizes: (" + << A.m() << ", " << A.n() << ") != (" + << B.m() << ", " << B.n() << ')' << nl << abort(FatalError); } Form AB(A.m(), A.n()); - Type* ABv = AB.v(); - const Type* Av = A.v(); - const Type* Bv = B.v(); + Type* ABv = AB.data(); + const Type* Av = A.cdata(); + const Type* Bv = B.cdata(); - const label mn = A.size(); - for (label i=0; i& A, const Matrix& B) template -Form Foam::operator*(const scalar s, const Matrix& M) +Form Foam::operator*(const scalar s, const Matrix& mat) { - Form sM(M.m(), M.n()); + Form result(mat.m(), mat.n()); - if (M.m() && M.n()) + const label len = mat.size(); + + if (len) { - Type* sMv = sM.v(); - const Type* Mv = M.v(); + Type* out = result.data(); + const Type* in = mat.cdata(); - const label mn = M.size(); - for (label i=0; i -Form Foam::operator*(const Matrix& M, const scalar s) +Form Foam::operator*(const Matrix& mat, const scalar s) { - Form sM(M.m(), M.n()); + Form result(mat.m(), mat.n()); - if (M.m() && M.n()) + const label len = mat.size(); + + if (len) { - Type* sMv = sM.v(); - const Type* Mv = M.v(); + Type* out = result.data(); + const Type* in = mat.cdata(); - const label mn = M.size(); - for (label i=0; i -Form Foam::operator/(const Matrix& M, const scalar s) +Form Foam::operator/(const Matrix& mat, const scalar s) { - Form sM(M.m(), M.n()); + Form result(mat.m(), mat.n()); - if (M.m() && M.n()) + const label len = mat.size(); + + if (len) { - Type* sMv = sM.v(); - const Type* Mv = M.v(); + Type* out = result.data(); + const Type* in = mat.cdata(); - const label mn = M.size(); - for (label i=0; i inline Foam::tmp> Foam::operator* ( - const Matrix& M, - const Field& f + const Matrix& mat, + const Field& x ) { - if (M.n() != f.size()) + if (mat.n() != x.size()) { FatalErrorInFunction << "Attempt to multiply incompatible matrix and field:" << nl - << "Matrix : " << M.m() << " x " << M.n() << nl - << "Field : " << f.size() << " rows" << nl - << "In order to multiply a matrix M and field f, " - "columns of M must equal rows of f" + << "Matrix : (" << mat.m() << ", " << mat.n() << ')' << nl + << "Field : " << x.size() << " rows" << nl + << "The number of matrix columns must equal the field size" << nl << abort(FatalError); } - tmp> tMf(new Field(M.m(), Zero)); - Field& Mf = tMf.ref(); + auto tresult = tmp>::New(mat.m(), Zero); + Field& result = tresult.ref(); - for (label i=0; i. + The layout is (mRows x nCols) - row-major order: + + \verbatim + (0,0) + +---> j [nCols] + | + | + v + i [mRows] + \endverbatim SourceFiles Matrix.C @@ -39,11 +49,8 @@ SourceFiles #ifndef Matrix_H #define Matrix_H -#include "bool.H" -#include "label.H" #include "uLabel.H" #include "Field.H" -#include "zero.H" #include "autoPtr.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -51,27 +58,10 @@ SourceFiles namespace Foam { -// Forward declaration of friend functions and operators - +// Forward Declarations template class Matrix; - -template Istream& operator>> -( - Istream&, - Matrix& -); - -template Ostream& operator<< -( - Ostream&, - const Matrix& -); - -template -class ConstMatrixBlock; - -template -class MatrixBlock; +template class ConstMatrixBlock; +template class MatrixBlock; /*---------------------------------------------------------------------------*\ @@ -81,16 +71,16 @@ class MatrixBlock; template class Matrix { - // Private data + // Private Data //- Number of rows and columns in Matrix. label mRows_, nCols_; - //- Row pointers + //- Vector of values of type Type Type* __restrict__ v_; - //- Allocate the storage for the element vector - void allocate(); + //- Allocate storage for the contents + inline void doAlloc(); public: @@ -105,7 +95,16 @@ public: // Static Member Functions //- Return a null Matrix - inline static const mType& null(); + inline static const Matrix& null(); + + + // Iterators + + //- Random access iterator for traversing a Matrix + typedef Type* iterator; + + //- Random access iterator for traversing a Matrix + typedef const Type* const_iterator; // Constructors @@ -117,30 +116,33 @@ public: Matrix(const label m, const label n); //- Construct with given number of rows and columns - // initializing all elements to zero + //- initializing all elements to zero Matrix(const label m, const label n, const zero); //- Construct with given number of rows and columns - // initializing all elements to the given value - Matrix(const label m, const label n, const Type&); + //- initializing all elements to the given value + Matrix(const label m, const label n, const Type& val); - //- Copy constructor. - Matrix(const mType&); + //- Copy construct + Matrix(const Matrix& mat); + + //- Move construct + Matrix(Matrix&& mat); //- Copy constructor from matrix of a different form template - explicit Matrix(const Matrix&); + explicit Matrix(const Matrix& mat); //- Construct from a block of another matrix template - Matrix(const ConstMatrixBlock&); + Matrix(const ConstMatrixBlock& Mb); //- Construct from a block of another matrix template - Matrix(const MatrixBlock&); + Matrix(const MatrixBlock& Mb); //- Construct from Istream. - Matrix(Istream&); + explicit Matrix(Istream& is); //- Clone inline autoPtr clone() const; @@ -152,25 +154,47 @@ public: // Member Functions - // Access + // Access - //- Return the number of rows - inline label m() const; + //- Return the number of rows. + inline label m() const noexcept; - //- Return the number of columns - inline label n() const; + //- Return the number of columns. + inline label n() const noexcept; - //- Return the number of elements in matrix (m*n) - inline label size() const; + //- Return the number of elements in matrix (m*n) + inline label size() const; - //- Return element vector of the constant Matrix - inline const Type* v() const; + //- Return true if the matrix is empty (ie, size() is zero) + inline bool empty() const noexcept; - //- Return element vector of the Matrix - inline Type* v(); + //- Return const pointer to the first data element, which can also + //- be used to address into the matrix contents + inline const Type* cdata() const noexcept; + + //- Return pointer to the first data element, which can also + //- be used to address into the matrix contents + inline Type* data() noexcept; + + //- Return const pointer to data in the specified row. + // Subscript checking only with FULLDEBUG + inline const Type* rowData(const label irow) const; + + //- Return pointer to data in the specified row. + // Subscript checking only with FULLDEBUG + inline Type* rowData(const label irow); + + //- Linear addressing const element access + // Subscript checking only with FULLDEBUG + inline const Type& at(const label idx) const; + + //- Linear addressing element access + // Subscript checking only with FULLDEBUG + inline Type& at(const label idx); - // Block access + + // Block Access inline ConstMatrixBlock block ( @@ -229,145 +253,246 @@ public: const label nStart ); + // Check - // Check + //- Check index i is within valid range (0 ... m-1). + inline void checki(const label irow) const; - //- Check index i is within valid range (0 ... m-1). - inline void checki(const label i) const; + //- Check index j is within valid range (0 ... n-1). + inline void checkj(const label jcol) const; - //- Check index j is within valid range (0 ... n-1). - inline void checkj(const label j) const; + //- Check that dimensions are positive, non-zero + inline void checkSize() const; + + //- True if all entries have identical values, and matrix is non-empty + inline bool uniform() const; - // Edit + // Edit - //- Clear the Matrix, i.e. set sizes to zero. - void clear(); + //- Clear the Matrix, i.e. set sizes to zero. + void clear(); - //- Transfer the contents of the argument Matrix into this Matrix - // and annul the argument Matrix. - void transfer(mType&); + //- Release storage management of the Matrix contents by transferring + //- management to a List + List release(); - //- Resize the matrix preserving the elements - void setSize(const label m, const label n); + //- Swap contents + void swap(Matrix& mat); - //- Resize the matrix without reallocating storage (unsafe) - inline void shallowResize(const label m, const label n); + //- Transfer the contents of the argument Matrix into this Matrix + //- and annul the argument Matrix. + void transfer(Matrix& mat); + + //- Change the matrix dimensions, preserving the elements + void resize(const label m, const label n); + + //- Change the matrix dimensions, preserving the elements + inline void setSize(const label m, const label n); + + //- Resize the matrix without reallocating storage (unsafe) + inline void shallowResize(const label m, const label n); //- Return the transpose of the matrix Form T() const; - // Member operators + // Member Operators - //- Return subscript-checked row of Matrix. - inline Type* operator[](const label); + //- Return const pointer to data in the specified row - rowData(). + // Subscript checking only with FULLDEBUG + inline const Type* operator[](const label irow) const; - //- Return subscript-checked row of constant Matrix. - inline const Type* operator[](const label) const; + //- Return pointer to data in the specified row - rowData(). + // Subscript checking only with FULLDEBUG + inline Type* operator[](const label irow); //- (i, j) const element access operator - inline const Type& operator()(const label i, const label j) const; + // Subscript checking only with FULLDEBUG + inline const Type& operator()(const label irow, const label jcol) const; //- (i, j) element access operator - inline Type& operator()(const label i, const label j); + // Subscript checking only with FULLDEBUG + inline Type& operator()(const label irow, const label jcol); - //- Assignment operator. Takes linear time. - void operator=(const mType&); + //- Copy assignment. Takes linear time. + void operator=(const Matrix& mat); + + //- Move assignment + void operator=(Matrix&& mat); //- Assignment to a block of another matrix template - void operator=(const ConstMatrixBlock&); + void operator=(const ConstMatrixBlock& Mb); //- Assignment to a block of another matrix template - void operator=(const MatrixBlock&); + void operator=(const MatrixBlock& Mb); //- Assignment of all elements to zero void operator=(const zero); //- Assignment of all elements to the given value - void operator=(const Type&); + void operator=(const Type& val); + + //- Matrix addition + void operator+=(const Matrix& other); + + //- Matrix subtraction + void operator-=(const Matrix& other); + + //- Matrix scalar multiplication + void operator*=(const scalar s); + + //- Matrix scalar division + void operator/=(const scalar s); - // IOstream operators + // Random access iterator (non-const) - //- Read Matrix from Istream, discarding contents of existing Matrix. - friend Istream& operator>> - ( - Istream&, - mType& - ); + //- Return an iterator to begin traversing a Matrix + inline iterator begin(); - //- Write Matrix to Ostream. - friend Ostream& operator<< - ( - Ostream&, - const mType& - ); + //- Return an iterator to end traversing a Matrix + inline iterator end(); + + + // Random access iterator (const) + + //- Return const_iterator to begin traversing a constant Matrix + inline const_iterator cbegin() const; + + //- Return const_iterator to end traversing a constant Matrix + inline const_iterator cend() const; + + //- Return const_iterator to begin traversing a constant Matrix + inline const_iterator begin() const; + + //- Return const_iterator to end traversing a constant Matrix + inline const_iterator end() const; + + + // IO + + //- Read Matrix from Istream, discarding existing contents. + bool readMatrix(Istream& is); + + //- Write Matrix, with line-breaks in ASCII when length exceeds + //- shortLen. + // Using '0' suppresses line-breaks entirely. + Ostream& writeMatrix(Ostream& os, const label shortLen=0) const; + + + // Housekeeping + + //- Deprecated(2019-04) raw data pointer, const access + // \deprecated(2019-04) - use cdata() method + const Type* FOAM_DEPRECATED_FOR(2019-04, "cdata() method") v() const + { + return v_; + } + + //- Deprecated(2019-04) raw data pointer, non-const access + // \deprecated(2019-04) - use data() method + Type* FOAM_DEPRECATED_FOR(2019-04, "data() method") v() + { + return v_; + } }; -// Global functions and operators +// IOstream Operators +//- Read Matrix from Istream, discarding contents of existing Matrix. template -const Type& max(const Matrix&); +Istream& operator>>(Istream& is, Matrix& mat); +//- Write Matrix to Ostream, as per Matrix::writeMatrix() with +//- default length, which is given by Detail::ListPolicy::short_length template -const Type& min(const Matrix&); +Ostream& operator<<(Ostream& os, const Matrix& mat); + +// Global Functions, Operators + +//- Find max value in the matrix template -Form operator-(const Matrix&); +const Type& max(const Matrix& mat); +//- Find min value in the matrix +template +const Type& min(const Matrix& mat); + +//- Find the min/max values of the matrix +template +MinMax minMax(const Matrix& mat); + +//- Matrix negation +template +Form operator-(const Matrix& mat); + +//- Matrix addition template Form operator+ ( - const Matrix&, - const Matrix& + const Matrix& A, + const Matrix& B ); + +//- Matrix subtraction template Form operator- ( - const Matrix&, - const Matrix& + const Matrix& A, + const Matrix& B ); + +//- Scalar multiplication of a matrix template Form operator* ( - const scalar, - const Matrix& + const scalar s, + const Matrix& mat ); + +//- Scalar multiplication of a matrix template Form operator* ( - const Matrix&, - const scalar + const Matrix& mat, + const scalar s ); + +//- Scalar division of a matrix template Form operator/ ( - const Matrix&, - const scalar + const Matrix& mat, + const scalar s ); + +//- Matrix-matrix multiplication template typename typeOfInnerProduct::type operator* ( - const Matrix& a, - const Matrix& b + const Matrix& A, + const Matrix& B ); + +//- Matrix-vector multiplication template tmp> operator* ( - const Matrix&, - const Field& + const Matrix& mat, + const Field& x ); diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H index bcc4a85675..6b6519cb80 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -27,6 +27,20 @@ License #include "MatrixBlock.H" +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +template +inline void Foam::Matrix::doAlloc() +{ + const label len = size(); + + if (len) + { + v_ = new Type[len]; + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template @@ -39,8 +53,8 @@ inline Foam::Matrix::Matrix() template -inline Foam::autoPtr> Foam::Matrix:: -clone() const +inline Foam::autoPtr> +Foam::Matrix::clone() const { return autoPtr>::New(*this); } @@ -56,14 +70,14 @@ inline const Foam::Matrix& Foam::Matrix::null() template -inline Foam::label Foam::Matrix::m() const +inline Foam::label Foam::Matrix::m() const noexcept { return mRows_; } template -inline Foam::label Foam::Matrix::n() const +inline Foam::label Foam::Matrix::n() const noexcept { return nCols_; } @@ -72,7 +86,14 @@ inline Foam::label Foam::Matrix::n() const template inline Foam::label Foam::Matrix::size() const { - return mRows_*nCols_; + return mRows_ * nCols_; +} + + +template +inline bool Foam::Matrix::empty() const noexcept +{ + return !mRows_ || !nCols_; } @@ -86,7 +107,7 @@ inline void Foam::Matrix::checki(const label i) const << "Attempt to access element from empty matrix" << abort(FatalError); } - else if (i<0 || i>=mRows_) + if (i < 0 || i >= mRows_) { FatalErrorInFunction << "Index " << i << " out of range 0 ... " << mRows_-1 @@ -106,7 +127,7 @@ inline void Foam::Matrix::checkj(const label j) const << "Attempt to access element from empty matrix" << abort(FatalError); } - else if (j<0 || j>=nCols_) + if (j < 0 || j >= nCols_) { FatalErrorInFunction << "index " << j << " out of range 0 ... " << nCols_-1 @@ -117,19 +138,104 @@ inline void Foam::Matrix::checkj(const label j) const template -inline const Type* Foam::Matrix::v() const +inline void Foam::Matrix::checkSize() const +{ + if (mRows_ < 0 || nCols_ < 0) + { + FatalErrorInFunction + << "Incorrect size (" << mRows_ << ", " << nCols_ << ')' << nl + << abort(FatalError); + } + // Could also check for odd sizes, like (0, N) and make (0,0) +} + + +template +inline bool Foam::Matrix::uniform() const +{ + const label len = size(); + + if (len == 0) + { + return false; + } + + for (label idx=1; idx +inline const Type* Foam::Matrix::cdata() const noexcept { return v_; } template -inline Type* Foam::Matrix::v() +inline Type* Foam::Matrix::data() noexcept { return v_; } +template +inline const Type* Foam::Matrix::rowData(const label irow) const +{ + #ifdef FULLDEBUG + checki(irow); + #endif + return (v_ + irow*nCols_); +} + + +template +inline Type* Foam::Matrix::rowData(const label irow) +{ + #ifdef FULLDEBUG + checki(irow); + #endif + return (v_ + irow*nCols_); +} + + +template +inline const Type& Foam::Matrix::at(const label idx) const +{ + #ifdef FULLDEBUG + if (idx < 0 || idx >= this->size()) + { + FatalErrorInFunction + << "index " << idx << " out of range 0 ... " << this->size() + << abort(FatalError); + } + #endif + return (v_ + idx); +} + + +template +inline Type& Foam::Matrix::at(const label idx) +{ + #ifdef FULLDEBUG + if (idx < 0 || idx >= this->size()) + { + FatalErrorInFunction + << "index " << idx << " out of range 0 ... " << this->size() + << abort(FatalError); + } + #endif + return (v_ + idx); +} + + template inline Foam::ConstMatrixBlock> Foam::Matrix::block @@ -282,6 +388,13 @@ Foam::Matrix::col } +template +inline void Foam::Matrix::setSize(const label m, const label n) +{ + resize(m, n); +} + + template void Foam::Matrix::shallowResize(const label m, const label n) { @@ -290,47 +403,97 @@ void Foam::Matrix::shallowResize(const label m, const label n) } +// * * * * * * * * * * * * * * * * * Iterators * * * * * * * * * * * * * * * // + +template +inline typename Foam::Matrix::iterator +Foam::Matrix::begin() +{ + return v_; +} + + +template +inline typename Foam::Matrix::iterator +Foam::Matrix::end() +{ + return v_ + (mRows_ * nCols_); +} + + +template +inline typename Foam::Matrix::const_iterator +Foam::Matrix::cbegin() const +{ + return v_; +} + + +template +inline typename Foam::Matrix::const_iterator +Foam::Matrix::cend() const +{ + return v_ + (mRows_ * nCols_); +} + + +template +inline typename Foam::Matrix::const_iterator +Foam::Matrix::begin() const +{ + return v_; +} + + +template +inline typename Foam::Matrix::const_iterator +Foam::Matrix::end() const +{ + return v_ + (mRows_ * nCols_); +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template inline const Type& Foam::Matrix::operator() ( - const label i, - const label j + const label irow, + const label jcol ) const { - checki(i); - checkj(j); - return v_[i*nCols_ + j]; + checki(irow); + checkj(jcol); + return v_[irow*nCols_ + jcol]; } template inline Type& Foam::Matrix::operator() ( - const label i, - const label j + const label irow, + const label jcol ) { - checki(i); - checkj(j); - return v_[i*nCols_ + j]; + checki(irow); + checkj(jcol); + return v_[irow*nCols_ + jcol]; } template -inline const Type* Foam::Matrix::operator[](const label i) const +inline const Type* Foam::Matrix::operator[](const label irow) const { - checki(i); - return v_ + i*nCols_; + checki(irow); + return v_ + irow*nCols_; } template -inline Type* Foam::Matrix::operator[](const label i) +inline Type* Foam::Matrix::operator[](const label irow) { - checki(i); - return v_ + i*nCols_; + checki(irow); + return v_ + irow*nCols_; } diff --git a/src/OpenFOAM/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C index e4628fbc09..67f5cfcc0b 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixIO.C +++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -30,6 +30,8 @@ License #include "Ostream.H" #include "token.H" #include "contiguous.H" +#include "ListPolicy.H" +#include // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -45,26 +47,25 @@ Foam::Matrix::Matrix(Istream& is) template -Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) +bool Foam::Matrix::readMatrix(Istream& is) { // Anull matrix - M.clear(); + clear(); is.fatalCheck(FUNCTION_NAME); token firstToken(is); - is.fatalCheck - ( - "operator>>(Istream&, Matrix&) : reading first token" - ); + is.fatalCheck("readMatrix : reading first token"); if (firstToken.isLabel()) { - M.mRows_ = firstToken.labelToken(); - M.nCols_ = readLabel(is); + mRows_ = firstToken.labelToken(); + nCols_ = readLabel(is); + doAlloc(); - label mn = M.mRows_*M.nCols_; + // The total size + const label len = size(); // Read list contents depending on data format if (is.format() == IOstream::ASCII || !contiguous()) @@ -72,29 +73,22 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) // Read beginning of contents char listDelimiter = is.readBeginList("Matrix"); - if (mn) + if (len) { - M.allocate(); - Type* v = M.v_; - if (listDelimiter == token::BEGIN_LIST) { - label k = 0; + label idx = 0; - // loop over rows - for (label i=0; i> v[k++]; + is >> v_[idx++]; - is.fatalCheck - ( - "operator>>(Istream&, Matrix&) : " - "reading entry" - ); + is.fatalCheck("readMatrix : reading reading entry"); } is.readEndList("MatrixRow"); @@ -105,16 +99,9 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) Type element; is >> element; - is.fatalCheck - ( - "operator>>(Istream&, Matrix&) : " - "reading the single entry" - ); + is.fatalCheck("readMatrix : reading the single entry"); - for (label i=0; i>(Istream& is, Matrix& M) } else { - if (mn) + if (len) { - M.allocate(); - Type* v = M.v_; + is.read(reinterpret_cast(v_), len*sizeof(Type)); - is.read(reinterpret_cast(v), mn*sizeof(Type)); - - is.fatalCheck - ( - "operator>>(Istream&, Matrix&) : " - "reading the binary block" - ); + is.fatalCheck("readMatrix : reading the binary block"); } } - } - else - { - FatalIOErrorInFunction(is) - << "incorrect first token, expected , found " - << firstToken.info() - << exit(FatalIOError); + + return len; } - return is; + + FatalIOErrorInFunction(is) + << "incorrect first token, expected , found " + << firstToken.info() + << exit(FatalIOError); + + return 0; } template -Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) +Foam::Ostream& Foam::Matrix::writeMatrix +( + Ostream& os, + const label shortLen +) const { - label mn = M.mRows_*M.nCols_; + const Matrix& mat = *this; - os << M.m() << token::SPACE << M.n(); + // The total size + const label len = mat.size(); + + // Rows, columns size + os << mat.m() << token::SPACE << mat.n(); // Write list contents depending on data format if (os.format() == IOstream::ASCII || !contiguous()) { - if (mn) + if (len) { - const Type* v = M.v_; + const Type* v = mat.cdata(); - // can the contents be considered 'uniform' (ie, identical) - bool uniform = (mn > 1 && contiguous()); - if (uniform) + // Can the contents be considered 'uniform' (ie, identical) + if (len > 1 && contiguous() && mat.uniform()) { - for (label i=0; i()) + else if (len < shortLen && contiguous()) { // Write start contents delimiter os << token::BEGIN_LIST; - label k = 0; + label idx = 0; - // loop over rows - for (label i=0; i& M) // Write start contents delimiter os << nl << token::BEGIN_LIST; - label k = 0; + label idx = 0; - // loop over rows - for (label i=0; i& M) } else { - if (mn) + // Contents are binary and contiguous + + if (len) { - os.write(reinterpret_cast(M.v_), mn*sizeof(Type)); + // write(...) includes surrounding start/end delimiters + os.write + ( + reinterpret_cast(mat.cdata()), + len*sizeof(Type) + ); } } @@ -257,4 +234,19 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) } +template +Foam::Istream& Foam::operator>>(Istream& is, Matrix& mat) +{ + mat.readMatrix(is); + return is; +} + + +template +Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& mat) +{ + return mat.writeMatrix(os, Detail::ListPolicy::short_length::value); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C index 75cd80a052..1d754b5a19 100644 --- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C +++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -30,22 +30,22 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::QRMatrix::QRMatrix(const MatrixType& M) +Foam::QRMatrix::QRMatrix(const MatrixType& mat) { - decompose(M); + decompose(mat); } // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // template -void Foam::QRMatrix::decompose(const MatrixType& M) +void Foam::QRMatrix::decompose(const MatrixType& mat) { - const label m = M.m(); - const label n = M.n(); + const label m = mat.m(); + const label n = mat.n(); // Initialize the R-matrix to M - R_ = M; + R_ = mat; // Initialize the Q-matrix to I Q_.setSize(m); @@ -158,9 +158,10 @@ void Foam::QRMatrix::decompose(const MatrixType& M) template +template class ListContainer> void Foam::QRMatrix::solvex ( - Field& x + ListContainer& x ) const { const label n = R_.n(); @@ -189,10 +190,11 @@ void Foam::QRMatrix::solvex template void Foam::QRMatrix::solve ( - Field& x, - const Field& source + List& x, + const UList& source ) const { + // Assert (&x != &source) ? const label m = Q_.m(); // x = Q_.T()*source; @@ -213,15 +215,14 @@ template Foam::tmp> Foam::QRMatrix::solve ( - const Field& source + const UList& source ) const { - tmp> tx(new Field(Q_.m())); - Field& x = tx.ref(); + auto tresult(tmp>::New(Q_.m())); - solve(x, source); + solve(tresult.ref(), source); - return tx; + return tresult; } diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H index e925967646..b716cd1954 100644 --- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H +++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -53,16 +53,16 @@ namespace Foam template class QRMatrix { - public: typedef typename MatrixType::cmptType cmptType; typedef SquareMatrix QMatrixType; typedef MatrixType RMatrixType; + private: - // Private data + // Private Data //- The Q-matrix QMatrixType Q_; @@ -71,12 +71,13 @@ private: RMatrixType R_; - // Private member functions + // Private Member Functions //- Solve the linear system with the Field argument x initialized to - // the appropriate transformed source (e.g. Q.T()*source) - // and return the solution in x - void solvex(Field& x) const; + //- the appropriate transformed source (e.g. Q.T()*source) + //- and return the solution in x + template class ListContainer> + void solvex(ListContainer& x) const; public: @@ -102,14 +103,14 @@ public: void decompose(const MatrixType& M); //- Solve the linear system with the given source - // and returning the solution in the Field argument x - void solve(Field& x, const Field& source) const; + //- and return the solution in the Field argument x + void solve(List& x, const UList& source) const; //- Solve the linear system with the given source - // returning the solution - tmp> solve(const Field& source) const; + //- and return the solution + tmp> solve(const UList& source) const; - //- Return the inverse of a square matrix + //- Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source QMatrixType inv() const; }; diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H index 256c22bac4..c1d6a0826b 100644 --- a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H +++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H index 1145776012..5b5e68e8b6 100644 --- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H @@ -77,27 +77,30 @@ public: inline RectangularMatrix(const MatrixBlock&); //- Construct with given number of rows and columns - // initializing all elements to zero + //- initializing all elements to zero inline RectangularMatrix(const label m, const label n, const zero); //- Construct with given number of rows and columns - // and value for all elements. - inline RectangularMatrix(const label m, const label n, const Type&); + //- and value for all elements. + inline RectangularMatrix(const label m, const label n, const Type& val); //- Construct as copy of a square matrix - inline RectangularMatrix(const SquareMatrix&); + inline RectangularMatrix(const SquareMatrix& mat); //- Construct from Istream. - inline RectangularMatrix(Istream&); + inline RectangularMatrix(Istream& is); //- Clone inline autoPtr> clone() const; - // Member operators + // Member Operators //- Assignment of all elements to zero void operator=(const zero); + + //- Assignment of all elements to the given value + void operator=(const Type& val); }; diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H index a584248b59..26c28ccfac 100644 --- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H @@ -84,20 +84,20 @@ inline Foam::RectangularMatrix::RectangularMatrix ( const label m, const label n, - const Type& t + const Type& val ) : - Matrix, Type>(m, n, t) + Matrix, Type>(m, n, val) {} template inline Foam::RectangularMatrix::RectangularMatrix ( - const SquareMatrix& SM + const SquareMatrix& mat ) : - Matrix, Type>(SM) + Matrix, Type>(mat) {} @@ -125,6 +125,13 @@ void Foam::RectangularMatrix::operator=(const zero) } +template +void Foam::RectangularMatrix::operator=(const Type& val) +{ + Matrix, Type>::operator=(val); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H index e662044924..44893b92a6 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -47,10 +47,8 @@ SourceFiles namespace Foam { -// Forward declaration of friend functions and operators - -template -class RectangularMatrix; +// Forward Declarations +template class RectangularMatrix; /*---------------------------------------------------------------------------*\ @@ -82,23 +80,23 @@ public: inline SquareMatrix(const MatrixBlock&); //- Construct given number of rows/columns - // initializing all elements to zero + //- initializing all elements to zero inline SquareMatrix(const label n, const zero); //- Construct given number of rows and columns (checked to be equal) - // initializing all elements to zero + //- initializing all elements to zero inline SquareMatrix(const label m, const label n, const zero); //- Construct given number of rows/columns - // Initializing to the identity matrix + //- initializing to the identity matrix inline SquareMatrix(const label n, const Identity); //- Construct with given number of rows and rows - // initializing all elements to the given value - inline SquareMatrix(const label n, const Type&); + //- initializing all elements to the given value + inline SquareMatrix(const label n, const Type& val); //- Construct as copy of a RectangularMatrix - // which is checked to be square + //- which is checked to be square inline explicit SquareMatrix(const RectangularMatrix&); //- Construct from Istream. @@ -110,13 +108,16 @@ public: // Member Functions - // Edit + // Edit - //- Resize the matrix preserving the elements - inline void setSize(const label m); + //- Resize the matrix preserving the elements + inline void resize(const label m); - //- Resize the matrix without reallocating storage (unsafe) - inline void shallowResize(const label m); + //- Resize the matrix preserving the elements + inline void setSize(const label m); + + //- Resize the matrix without reallocating storage (unsafe) + inline void shallowResize(const label m); // Member operators @@ -124,7 +125,10 @@ public: //- Assignment of all elements to zero void operator=(const zero); - //- Assignment elements to the + //- Assignment of all elements to the given value + void operator=(const Type& val); + + //- Clear and assign diagonal to identity void operator=(const Identity); }; diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H index 529ed82ba8..0d5a3a8329 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -94,6 +94,17 @@ inline Foam::SquareMatrix::SquareMatrix } +template +inline Foam::SquareMatrix::SquareMatrix +( + const label n, + const Type& val +) +: + Matrix, Type>(n, n, val) +{} + + template inline Foam::SquareMatrix::SquareMatrix ( @@ -113,27 +124,16 @@ inline Foam::SquareMatrix::SquareMatrix template inline Foam::SquareMatrix::SquareMatrix ( - const label n, - const Type& t + const RectangularMatrix& mat ) : - Matrix, Type>(n, n, t) -{} - - -template -inline Foam::SquareMatrix::SquareMatrix -( - const RectangularMatrix& RM -) -: - Matrix, Type>(RM) + Matrix, Type>(mat) { - if (this->m() != this->n()) + if (mat.m() != mat.n()) { FatalErrorInFunction << "Attempt to construct a square matrix from a rectangular matrix " - << this->m() << " x " << this->n() << nl + << mat.m() << " x " << mat.n() << nl << abort(FatalError); } } @@ -156,10 +156,17 @@ Foam::SquareMatrix::clone() const // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +inline void Foam::SquareMatrix::resize(const label m) +{ + Matrix, Type>::resize(m, m); +} + + template inline void Foam::SquareMatrix::setSize(const label m) { - Matrix, Type>::setSize(m, m); + Matrix, Type>::resize(m, m); } @@ -179,10 +186,18 @@ void Foam::SquareMatrix::operator=(const zero) } +template +void Foam::SquareMatrix::operator=(const Type& val) +{ + Matrix, Type>::operator=(val); +} + + template void Foam::SquareMatrix::operator=(const Identity) { Matrix, Type>::operator=(Zero); + for (label i=0; in(); i++) { this->operator()(i, i) = Type(I); diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H index 267da6b327..f4bf7c9c41 100644 --- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H +++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -70,13 +70,13 @@ public: //- Construct given number of rows/columns, initializing to zero inline SymmetricSquareMatrix(const label n, const zero); + //- Construct with given number of rows/columns + //- initializing all elements to the given value + inline SymmetricSquareMatrix(const label n, const Type& val); + //- Construct given number of rows/columns, inline SymmetricSquareMatrix(const label n, const Identity); - //- Construct with given number of rows/columns - // initializing all elements to the given value - inline SymmetricSquareMatrix(const label n, const Type&); - //- Construct from Istream. inline SymmetricSquareMatrix(Istream&); @@ -85,7 +85,7 @@ public: }; -// Global functions +// Global Functions //- Return the LU decomposed SymmetricSquareMatrix inverse template diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H index 7dcc04e2b7..a2f8dc0109 100644 --- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H +++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -52,6 +52,17 @@ inline Foam::SymmetricSquareMatrix::SymmetricSquareMatrix {} +template +inline Foam::SymmetricSquareMatrix::SymmetricSquareMatrix +( + const label n, + const Type& val +) +: + Matrix, Type>(n, n, val) +{} + + template inline Foam::SymmetricSquareMatrix::SymmetricSquareMatrix ( @@ -68,17 +79,6 @@ inline Foam::SymmetricSquareMatrix::SymmetricSquareMatrix } -template -inline Foam::SymmetricSquareMatrix::SymmetricSquareMatrix -( - const label n, - const Type& t -) -: - Matrix, Type>(n, n, t) -{} - - template inline Foam::SymmetricSquareMatrix::SymmetricSquareMatrix(Istream& is) :