diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index bf70af621fb6a4a466929bf364a00338eb35aab3..61aaf2d8a35ae0c89c3f9f38a8e7004cf4b183a7 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 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<scalar> hmm(3);
+ // SquareMatrix
+ {
+ Info<< nl << "Test SquareMatrix" << nl;
+
+ SquareMatrix<scalar> 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<scalar> stole(square1.release());
+
+ Info<< "matrix: " << square1
+ << " begin: " << long(square1.cdata()) << nl;
+
+ Info<< "List: " << stole << nl;
- 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<< "min/max: " << minMax(square1) << nl;
- //Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
- Info<< max(hmm) << endl;
- Info<< min(hmm) << endl;
+ square1 = 100;
- SquareMatrix<scalar> hmm2(3, I);
+ Info<< "matrix: " << square1
+ << " begin: " << long(square1.cdata()) << nl;
- hmm = hmm2;
- Info<< hmm << endl;
+ SquareMatrix<scalar> square2(3, I);
- //SquareMatrix<scalar> hmm3(Sin);
+ square1 = square2;
- //Info<< hmm3 << endl;
+ Info<< nl << "After copy assign from Identity:" << nl << square1 << nl;
- SquareMatrix<scalar> hmm4;
+ square1 += 1.25*square2;
- hmm4 = hmm2;
+ Info<< nl << "After +=" << nl << square1 << nl;
- Info<< hmm4 << endl;
+ square1 -= square2.T();
- SquareMatrix<scalar> hmm5;
+ Info<< nl << "After -=" << nl << square1 << nl;
- hmm4 = hmm5;
- Info<< hmm5 << endl;
+ square1 *= 10;
+ Info<< nl << "After *=" << nl << square1 << nl;
+
+ square1 /= 8;
+
+ Info<< nl << "After /=" << nl << square1 << nl;
+
+ SquareMatrix<scalar> square4;
+ square4 = square2;
+
+ Info<< nl << square4 << endl;
+
+ SquareMatrix<scalar> square5;
+ square4 = square5;
+ Info<< nl << square5 << endl;
+ }
+
+ // RectangularMatrix
{
+ Info<< nl << "Test RectangularMatrix" << nl;
+
RectangularMatrix<scalar> rm1(5, 6, 3.1);
rm1(0, 1) = 4.5;
RectangularMatrix<scalar> rm1b(rm1.block(2, 2, 0, 0));
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
index 1a3da1e7e7bfe75d1069e76b70accab35ace156b..99b7b1135993e81cbb19d34734a593633b2b5ec9 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -37,45 +37,54 @@ inline Foam::DiagonalMatrix<Type>::DiagonalMatrix()
template<class Type>
-template<class Form>
-Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n)
:
- List<Type>(min(a.m(), a.n()))
-{
- forAll(*this, i)
- {
- this->operator[](i) = a(i, i);
- }
-}
+ List<Type>(n)
+{}
template<class Type>
-Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size)
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n, const zero)
:
- List<Type>(size)
+ List<Type>(n, Zero)
{}
template<class Type>
-Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size, const Type& val)
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n, const Type& val)
:
- List<Type>(size, val)
+ List<Type>(n, val)
{}
+template<class Type>
+template<class Form>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& mat)
+:
+ List<Type>(min(mat.m(), mat.n()))
+{
+ label i = 0;
+
+ for (Type& val : *this)
+ {
+ val = mat(i, i);
+ ++i;
+ }
+}
+
+
template<class Type>
Foam::DiagonalMatrix<Type>& Foam::DiagonalMatrix<Type>::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<Type>& Foam::DiagonalMatrix<Type>::invert()
template<class Type>
-Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& A)
+Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& mat)
{
- DiagonalMatrix<Type> Ainv = A;
+ DiagonalMatrix<Type> 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 e0474803cb6254763779ccf5e025c2c9c5d82490..c5e56dac95e695182dfb35f307a04e9f736d4fd7 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 | Copyright (C) 2004-2011 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2011, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -27,8 +27,7 @@ Class
Foam::DiagonalMatrix
Description
- DiagonalMatrix<Type> is a 2D diagonal matrix of objects
- of type Type, size nxn
+ A 2D diagonal matrix of objects of type \<Type\>, size (N x N)
SourceFiles
DiagonalMatrix.C
@@ -45,12 +44,11 @@ SourceFiles
namespace Foam
{
-// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * //
-
+// Forward Declarations
template<class Form, class Type> class Matrix;
/*---------------------------------------------------------------------------*\
- Class DiagonalMatrix Declaration
+ Class DiagonalMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
@@ -62,24 +60,29 @@ public:
// Constructors
- //- Null constructor.
+ //- Construct null.
DiagonalMatrix<Type>();
- //- Construct from diagonal component of a Matrix
- template<class Form>
- DiagonalMatrix<Type>(const Matrix<Form, Type>&);
-
//- Construct empty from size
- DiagonalMatrix<Type>(const label size);
+ explicit DiagonalMatrix<Type>(const label n);
+
+ //- Construct from size
+ //- initializing all elements to the given value
+ DiagonalMatrix<Type>(const label n, const zero);
//- Construct from size and a value
- DiagonalMatrix<Type>(const label, const Type&);
+ DiagonalMatrix<Type>(const label n, const Type& val);
+
+ //- Construct from the diagonal of a Matrix
+ template<class Form>
+ DiagonalMatrix<Type>(const Matrix<Form, Type>& mat);
- // Member functions
+ // Member Functions
//- Invert the diagonal matrix and return itself
DiagonalMatrix<Type>& invert();
+
};
@@ -87,7 +90,7 @@ public:
//- Return the diagonal Matrix inverse
template<class Type>
-DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>&);
+DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>& mat);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
index c18bbe5c519179d51dc523c6ef64536f0d02084d..7798d2ea9dd7d9441be65ea1d998aaab61f0a6eb 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<Type>::LLTMatrix()
template<class Type>
-Foam::LLTMatrix<Type>::LLTMatrix(const SquareMatrix<Type>& M)
+Foam::LLTMatrix<Type>::LLTMatrix(const SquareMatrix<Type>& mat)
{
- decompose(M);
+ decompose(mat);
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class Type>
-void Foam::LLTMatrix<Type>::decompose(const SquareMatrix<Type>& M)
+void Foam::LLTMatrix<Type>::decompose(const SquareMatrix<Type>& mat)
{
SquareMatrix<Type>& LLT = *this;
// Initialize the LLT decomposition matrix to M
- LLT = M;
+ LLT = mat;
const label m = LLT.m();
- for (label i=0; i<m; i++)
+ for (label i=0; i<m; ++i)
{
- for (label j=0; j<m; j++)
+ for (label j=0; j<m; ++j)
{
if (j > i)
{
@@ -65,7 +65,7 @@ void Foam::LLTMatrix<Type>::decompose(const SquareMatrix<Type>& M)
Type sum = LLT(i, j);
- for (label k=0; k<j; k++)
+ for (label k=0; k<j; ++k)
{
sum -= LLT(i, k)*LLT(j, k);
}
@@ -93,11 +93,11 @@ void Foam::LLTMatrix<Type>::decompose(const SquareMatrix<Type>& M)
template<class Type>
void Foam::LLTMatrix<Type>::solve
(
- Field<Type>& x,
- const Field<Type>& source
+ List<Type>& x,
+ const UList<Type>& 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<Type>::solve
const SquareMatrix<Type>& LLT = *this;
const label m = LLT.m();
- for (label i=0; i<m; i++)
+ for (label i=0; i<m; ++i)
{
Type sum = source[i];
- for (label j=0; j<i; j++)
+ for (label j=0; j<i; ++j)
{
sum = sum - LLT(i, j)*x[j];
}
@@ -118,33 +118,31 @@ void Foam::LLTMatrix<Type>::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<m; j++)
+ for (label j=i + 1; j<m; ++j)
{
sum = sum - LLT(j, i)*x[j];
}
x[i] = sum/LLT(i, i);
}
-
}
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::LLTMatrix<Type>::solve
(
- const Field<Type>& source
+ const UList<Type>& source
) const
{
- tmp<Field<Type>> tx(new Field<Type>(this->m()));
- Field<Type>& x = tx.ref();
+ auto tresult(tmp<Field<Type>>::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 80d42c23b22e20143bdca8b36fd2ae8adc7c3ebc..a9f853e41ce43f2884389b18deaab35d21e8b592 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<Type>
{
-
public:
// Constructors
@@ -66,23 +65,23 @@ public:
//- Construct null
LLTMatrix();
- //- Construct from a square matrix and perform the decomposition
- LLTMatrix(const SquareMatrix<Type>& M);
+ //- Construct and perform the decomposition on input square matrix
+ LLTMatrix(const SquareMatrix<Type>& mat);
// Member Functions
- //- Perform the Cholesky decomposition of the matrix
- void decompose(const SquareMatrix<Type>& M);
+ //- Copy matrix and perform Cholesky decomposition
+ void decompose(const SquareMatrix<Type>& 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<Type>& x, const Field<Type>& source) const;
+ void solve(List<Type>& x, const UList<Type>& source) const;
//- Solve the linear system with the given source
- // returning the solution
- tmp<Field<Type>> solve(const Field<Type>& source) const;
+ //- returning the solution
+ tmp<Field<Type>> solve(const UList<Type>& source) const;
};
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
index 6a666f099fc7991253d544523cb8bb878e9ab7bc..c5fb4c48e46d265cd59481e10ad129dec70ab509 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<class Type>
- void solve(Field<Type>& x, const Field<Type>& source) const;
+ void solve(List<Type>& x, const UList<Type>& source) const;
//- Solve the linear system with the given source
// returning the solution
template<class Type>
- tmp<Field<Type>> solve(const Field<Type>& source) const;
+ tmp<Field<Type>> solve(const UList<Type>& 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 dc5b508d3d919e7b79a8b36585f1830a6fe9f2f0..30124aa7fec3be4ad8b0f39c5688326ff5515dba 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 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<class Type>
void Foam::LUscalarMatrix::solve
(
- Field<Type>& x,
- const Field<Type>& source
+ List<Type>& x,
+ const UList<Type>& source
) const
{
// If x and source are different initialize x = source
@@ -45,15 +45,13 @@ void Foam::LUscalarMatrix::solve
if (Pstream::parRun())
{
- Field<Type> X(m());
+ List<Type> X; // scratch space (on master)
if (Pstream::master(comm_))
{
- typename Field<Type>::subField
- (
- X,
- x.size()
- ) = x;
+ X.resize(m());
+
+ SubList<Type>(X, x.size()) = x;
for
(
@@ -82,7 +80,7 @@ void Foam::LUscalarMatrix::solve
(
Pstream::commsTypes::scheduled,
Pstream::masterNo(),
- reinterpret_cast<const char*>(x.begin()),
+ reinterpret_cast<const char*>(x.cdata()),
x.byteSize(),
Pstream::msgType(),
comm_
@@ -93,11 +91,7 @@ void Foam::LUscalarMatrix::solve
{
LUBacksubstitute(*this, pivotIndices_, X);
- x = typename Field<Type>::subField
- (
- X,
- x.size()
- );
+ x = SubList<Type>(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<char*>(x.begin()),
+ reinterpret_cast<char*>(x.data()),
x.byteSize(),
Pstream::msgType(),
comm_
@@ -143,13 +137,12 @@ void Foam::LUscalarMatrix::solve
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::LUscalarMatrix::solve
(
- const Field<Type>& source
+ const UList<Type>& source
) const
{
- tmp<Field<Type>> tx(new Field<Type>(m()));
- Field<Type>& x = tx.ref();
+ auto tx(tmp<Field<Type>>::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 d5d205a9a0dce1b850d0823e9b03f30f8e664912..4b3c44b8ebf1e8d091b766c38f71701435a0fcc9 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -26,18 +26,8 @@ License
\*---------------------------------------------------------------------------*/
#include "Matrix.H"
-
-// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::allocate()
-{
- if (mRows_ && nCols_)
- {
- v_ = new Type[size()];
- }
-}
-
+#include <functional>
+#include <algorithm>
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -48,14 +38,9 @@ Foam::Matrix<Form, Type>::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<Form, Type>::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<mn; i++)
- {
- v_[i] = Zero;
- }
- }
+ std::fill(begin(), end(), Zero);
}
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& s)
+Foam::Matrix<Form, Type>::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_)
- {
- const label mn = size();
- for (label i=0; i<mn; i++)
- {
- v_[i] = s;
- }
- }
+ std::fill(begin(), end(), val);
}
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& M)
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& mat)
:
- mRows_(M.mRows_),
- nCols_(M.nCols_),
+ mRows_(mat.mRows_),
+ nCols_(mat.nCols_),
v_(nullptr)
{
- if (M.v_)
+ if (mat.cdata())
{
- allocate();
+ doAlloc();
- const label mn = size();
- for (label i=0; i<mn; i++)
- {
- v_[i] = M.v_[i];
- }
+ std::copy(mat.cbegin(), mat.cend(), v_);
}
}
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(Matrix<Form, Type>&& mat)
+:
+ mRows_(mat.mRows_),
+ nCols_(mat.nCols_),
+ v_(mat.v_)
+{
+ mat.mRows_ = 0;
+ mat.nCols_ = 0;
+ mat.v_ = nullptr;
+}
+
+
template<class Form, class Type>
template<class Form2>
-Foam::Matrix<Form, Type>::Matrix(const Matrix<Form2, Type>& M)
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form2, Type>& 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<mn; i++)
- {
- v_[i] = M.v()[i];
- }
+ std::copy(mat.cbegin(), mat.cend(), v_);
}
}
@@ -164,13 +130,13 @@ inline Foam::Matrix<Form, Type>::Matrix
mRows_(Mb.m()),
nCols_(Mb.n())
{
- allocate();
+ doAlloc();
- for (label i=0; i<mRows_; i++)
+ for (label i=0; i < mRows_; ++i)
{
- for (label j=0; j<nCols_; j++)
+ for (label j=0; j < nCols_; ++j)
{
- (*this)(i,j) = Mb(i,j);
+ (*this)(i, j) = Mb(i,j);
}
}
}
@@ -186,13 +152,13 @@ inline Foam::Matrix<Form, Type>::Matrix
mRows_(Mb.m()),
nCols_(Mb.n())
{
- allocate();
+ doAlloc();
- for (label i=0; i<mRows_; i++)
+ for (label i=0; i < mRows_; ++i)
{
- for (label j=0; j<nCols_; j++)
+ for (label j=0; j<nCols_; ++j)
{
- (*this)(i,j) = Mb(i,j);
+ (*this)(i, j) = Mb(i, j);
}
}
}
@@ -227,32 +193,65 @@ void Foam::Matrix<Form, Type>::clear()
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& M)
+Foam::List<Type> Foam::Matrix<Form, Type>::release()
{
+ List<Type> list;
+
+ const label len = size();
+
+ if (v_ && len)
+ {
+ UList<Type> storage(v_, len);
+ list.swap(storage);
+
+ v_ = nullptr;
+ }
clear();
- mRows_ = M.mRows_;
- M.mRows_ = 0;
+ return list;
+}
- nCols_ = M.nCols_;
- M.nCols_ = 0;
- v_ = M.v_;
- M.v_ = nullptr;
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::swap(Matrix<Form, Type>& mat)
+{
+ Foam::Swap(mRows_, mat.mRows_);
+ Foam::Swap(nCols_, mat.nCols_);
+ Foam::Swap(v_, mat.v_);
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
+void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& mat)
{
- mType newMatrix(m, n, Zero);
+ clear();
+
+ mRows_ = mat.mRows_;
+ nCols_ = mat.nCols_;
+ v_ = mat.v_;
+
+ mat.mRows_ = 0;
+ mat.nCols_ = 0;
+ mat.v_ = nullptr;
+}
- label minM = min(m, mRows_);
- label minN = min(n, nCols_);
- for (label i=0; i<minM; i++)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::resize(const label m, const label n)
+{
+ if (m == mRows_ && n == nCols_)
+ {
+ return;
+ }
+
+ Matrix<Form, Type> 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<minN; j++)
+ for (label j=0; j < ncol; ++j)
{
newMatrix(i, j) = (*this)(i, j);
}
@@ -265,14 +264,13 @@ void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
template<class Form, class Type>
Form Foam::Matrix<Form, Type>::T() const
{
- const Matrix<Form, Type>& A = *this;
Form At(n(), m());
- for (label i=0; i<m(); i++)
+ for (label i=0; i < m(); ++i)
{
- for (label j=0; j<n(); j++)
+ for (label j=0; j < n(); ++j)
{
- At(j, i) = A(i, j);
+ At(j, i) = (*this)(i, j);
}
}
@@ -283,31 +281,41 @@ Form Foam::Matrix<Form, Type>::T() const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& M)
+void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& 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<mn; i++)
- {
- v_[i] = M.v_[i];
- }
+ std::copy(mat.cbegin(), mat.cend(), v_);
+ }
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(Matrix<Form, Type>&& mat)
+{
+ if (this == &mat)
+ {
+ FatalErrorInFunction
+ << "Attempted assignment to self"
+ << abort(FatalError);
}
+
+ this->transfer(mat);
}
@@ -318,11 +326,11 @@ void Foam::Matrix<Form, Type>::operator=
const ConstMatrixBlock<MatrixType>& Mb
)
{
- for (label i=0; i<mRows_; i++)
+ for (label i=0; i < mRows_; ++i)
{
- for (label j=0; j<nCols_; j++)
+ for (label j=0; j < nCols_; ++j)
{
- (*this)(i,j) = Mb(i,j);
+ (*this)(i, j) = Mb(i, j);
}
}
}
@@ -335,130 +343,170 @@ void Foam::Matrix<Form, Type>::operator=
const MatrixBlock<MatrixType>& Mb
)
{
- for (label i=0; i<mRows_; i++)
+ for (label i=0; i < mRows_; ++i)
{
- for (label j=0; j<nCols_; j++)
+ for (label j=0; j < nCols_; ++j)
{
- (*this)(i,j) = Mb(i,j);
+ (*this)(i, j) = Mb(i, j);
}
}
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Type& s)
+void Foam::Matrix<Form, Type>::operator=(const Type& val)
{
- if (v_)
- {
- const label mn = size();
- for (label i=0; i<mn; i++)
- {
- v_[i] = s;
- }
- }
+ std::fill(begin(), end(), val);
}
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator=(const zero)
{
- if (v_)
- {
- const label mn = size();
- for (label i=0; i<mn; i++)
- {
- v_[i] = Zero;
- }
- }
+ std::fill(begin(), end(), Zero);
}
-// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
-
template<class Form, class Type>
-const Type& Foam::max(const Matrix<Form, Type>& M)
+void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& other)
{
- const label mn = M.size();
+ if (this == &other)
+ {
+ FatalErrorInFunction
+ << "Attempted addition to self"
+ << abort(FatalError);
+ }
- if (mn)
+ if (m() != other.m() || n() != other.n())
{
- label curMaxI = 0;
- const Type* Mv = M.v();
+ FatalErrorInFunction
+ << "Attempt to add 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];
+ }
+}
- for (label i=1; i<mn; i++)
- {
- if (Mv[i] > Mv[curMaxI])
- {
- curMaxI = i;
- }
- }
- return Mv[curMaxI];
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator-=(const Matrix<Form, Type>& other)
+{
+ if (this == &other)
+ {
+ FatalErrorInFunction
+ << "Attempted subtraction from self"
+ << abort(FatalError);
}
- else
+
+ if (m() != other.m() || n() != other.n())
{
FatalErrorInFunction
- << "Matrix is empty"
+ << "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();
- // Return in error to keep compiler happy
- return M(0, 0);
+ for (label idx=0; idx < len; ++idx)
+ {
+ out[idx] -= in[idx];
}
}
template<class Form, class Type>
-const Type& Foam::min(const Matrix<Form, Type>& M)
+void Foam::Matrix<Form, Type>::operator*=(const scalar s)
{
- const label mn = M.size();
+ for (Type& val : *this)
+ {
+ val *= s;
+ }
+}
+
- if (mn)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator/=(const scalar s)
+{
+ for (Type& val : *this)
{
- label curMinI = 0;
- const Type* Mv = M.v();
+ val /= s;
+ }
+}
- for (label i=1; i<mn; i++)
- {
- if (Mv[i] < Mv[curMinI])
- {
- curMinI = i;
- }
- }
- return Mv[curMinI];
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
+
+template<class Form, class Type>
+const Type& Foam::max(const Matrix<Form, Type>& mat)
+{
+ if (mat.empty())
+ {
+ FatalErrorInFunction
+ << "Matrix is empty" << abort(FatalError);
}
- else
+
+ return *(std::max_element(mat.cbegin(), mat.cend()));
+}
+
+
+template<class Form, class Type>
+const Type& Foam::min(const Matrix<Form, Type>& mat)
+{
+ if (mat.empty())
{
FatalErrorInFunction
- << "Matrix is empty"
- << abort(FatalError);
+ << "Matrix is empty" << abort(FatalError);
+ }
+
+ return *(std::min_element(mat.cbegin(), mat.cend()));
+}
+
+
+template<class Form, class Type>
+Foam::MinMax<Type> Foam::minMax(const Matrix<Form, Type>& mat)
+{
+ MinMax<Type> result;
- // Return in error to keep compiler happy
- return M(0, 0);
+ for (const Type& val : mat)
+ {
+ result += val;
}
+
+ return result;
}
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& M)
+Form Foam::operator-(const Matrix<Form, Type>& mat)
{
- Form nM(M.m(), M.n());
-
- if (M.m() && M.n())
- {
- Type* nMv = nM.v();
- const Type* Mv = M.v();
+ Form result(mat.m(), mat.n());
- const label mn = M.size();
- for (label i=0; i<mn; i++)
- {
- nMv[i] = -Mv[i];
- }
- }
+ std::transform
+ (
+ mat.cbegin(),
+ mat.cend(),
+ result.begin(),
+ std::negate<Type>()
+ );
- return nM;
+ return result;
}
@@ -468,29 +516,23 @@ Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& 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<mn; i++)
+ const label len = A.size();
+
+ for (label idx=0; idx < len; ++idx)
{
- ABv[i] = Av[i] + Bv[i];
+ ABv[idx] = Av[idx] + Bv[idx];
}
return AB;
@@ -503,29 +545,23 @@ Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& 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 len = A.size();
- const label mn = A.size();
- for (label i=0; i<mn; i++)
+ for (label idx=0; idx < len; ++idx)
{
- ABv[i] = Av[i] - Bv[i];
+ ABv[idx] = Av[idx] - Bv[idx];
}
return AB;
@@ -533,65 +569,68 @@ Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
template<class Form, class Type>
-Form Foam::operator*(const scalar s, const Matrix<Form, Type>& M)
+Form Foam::operator*(const scalar s, const Matrix<Form, Type>& mat)
{
- Form sM(M.m(), M.n());
+ Form result(mat.m(), mat.n());
+
+ const label len = mat.size();
- if (M.m() && M.n())
+ 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<mn; i++)
+ for (label idx=0; idx < len; ++idx)
{
- sMv[i] = s*Mv[i];
+ out[idx] = s * in[idx];
}
}
- return sM;
+ return result;
}
template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& M, const scalar s)
+Form Foam::operator*(const Matrix<Form, Type>& 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<mn; i++)
+ for (label idx=0; idx < len; ++idx)
{
- sMv[i] = Mv[i]*s;
+ out[idx] = in[idx] * s;
}
}
- return sM;
+ return result;
}
template<class Form, class Type>
-Form Foam::operator/(const Matrix<Form, Type>& M, const scalar s)
+Form Foam::operator/(const Matrix<Form, Type>& mat, const scalar s)
{
- Form sM(M.m(), M.n());
+ Form result(mat.m(), mat.n());
+
+ const label len = mat.size();
- if (M.m() && M.n())
+ 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<mn; i++)
+ for (label idx=0; idx < len; ++idx)
{
- sMv[i] = Mv[i]/s;
+ out[idx] = in[idx] / s;
}
}
- return sM;
+ return result;
}
@@ -607,10 +646,9 @@ Foam::operator*
{
FatalErrorInFunction
<< "Attempt to multiply incompatible matrices:" << nl
- << "Matrix A : " << A.m() << " x " << A.n() << nl
- << "Matrix B : " << B.m() << " x " << B.n() << nl
- << "In order to multiply matrices, columns of A must equal "
- << "rows of B"
+ << "Matrix A : (" << A.m() << ", " << A.n() << ')' << nl
+ << "Matrix B : (" << B.m() << ", " << B.n() << ')' << nl
+ << "The columns of A must equal rows of B"
<< abort(FatalError);
}
@@ -621,13 +659,13 @@ Foam::operator*
Zero
);
- for (label i=0; i<AB.m(); i++)
+ for (label i=0; i < AB.m(); ++i)
{
- for (label j=0; j<AB.n(); j++)
+ for (label j=0; j < AB.n(); ++j)
{
- for (label k=0; k<B.m(); k++)
+ for (label k=0; k < B.m(); ++k)
{
- AB(i, j) += A(i, k)*B(k, j);
+ AB(i, j) += A(i, k) * B(k, j);
}
}
}
@@ -639,33 +677,32 @@ Foam::operator*
template<class Form, class Type>
inline Foam::tmp<Foam::Field<Type>> Foam::operator*
(
- const Matrix<Form, Type>& M,
- const Field<Type>& f
+ const Matrix<Form, Type>& mat,
+ const Field<Type>& 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<Field<Type>> tMf(new Field<Type>(M.m(), Zero));
- Field<Type>& Mf = tMf.ref();
+ auto tresult = tmp<Field<Type>>::New(mat.m(), Zero);
+ Field<Type>& result = tresult.ref();
- for (label i=0; i<M.m(); i++)
+ for (label i=0; i < mat.m(); ++i)
{
- for (label j=0; j<M.n(); j++)
+ for (label j=0; j < mat.n(); ++j)
{
- Mf[i] += M(i, j)*f[j];
+ result[i] += mat(i, j) * x[j];
}
}
- return tMf;
+ return tresult;
}
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 7c67a4b984c616d23588131fbebc479c80304fc9..700a31e523d303b805b2f1b30afdf08238e494d6 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -28,6 +28,16 @@ Class
Description
A templated (m x n) matrix of objects of \<T\>.
+ 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 Form, class Type> class Matrix;
-
-template<class Form, class Type> Istream& operator>>
-(
- Istream&,
- Matrix<Form, Type>&
-);
-
-template<class Form, class Type> Ostream& operator<<
-(
- Ostream&,
- const Matrix<Form, Type>&
-);
-
-template<class MatrixType>
-class ConstMatrixBlock;
-
-template<class MatrixType>
-class MatrixBlock;
+template<class MatrixType> class ConstMatrixBlock;
+template<class MatrixType> class MatrixBlock;
/*---------------------------------------------------------------------------*\
@@ -81,16 +71,16 @@ class MatrixBlock;
template<class Form, class Type>
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<Form, Type>& 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 construct
+ Matrix(const Matrix<Form, Type>& mat);
- //- Copy constructor.
- Matrix(const mType&);
+ //- Move construct
+ Matrix(Matrix<Form, Type>&& mat);
//- Copy constructor from matrix of a different form
template<class Form2>
- explicit Matrix(const Matrix<Form2, Type>&);
+ explicit Matrix(const Matrix<Form2, Type>& mat);
//- Construct from a block of another matrix
template<class MatrixType>
- Matrix(const ConstMatrixBlock<MatrixType>&);
+ Matrix(const ConstMatrixBlock<MatrixType>& Mb);
//- Construct from a block of another matrix
template<class MatrixType>
- Matrix(const MatrixBlock<MatrixType>&);
+ Matrix(const MatrixBlock<MatrixType>& Mb);
//- Construct from Istream.
- Matrix(Istream&);
+ explicit Matrix(Istream& is);
//- Clone
inline autoPtr<mType> 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;
- // Block access
+ //- 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
inline ConstMatrixBlock<mType> block
(
@@ -229,145 +253,246 @@ public:
const label nStart
);
+ // Check
+
+ //- Check index i is within valid range (0 ... m-1).
+ inline void checki(const label irow) const;
- // Check
+ //- Check index j is within valid range (0 ... n-1).
+ inline void checkj(const label jcol) const;
- //- Check index i is within valid range (0 ... m-1).
- inline void checki(const label i) const;
+ //- Check that dimensions are positive, non-zero
+ inline void checkSize() const;
- //- Check index j is within valid range (0 ... n-1).
- inline void checkj(const label j) 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<Type> release();
- //- Resize the matrix preserving the elements
- void setSize(const label m, const label n);
+ //- Swap contents
+ void swap(Matrix<Form, Type>& 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<Form, Type>& 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);
+
+ //- Copy assignment. Takes linear time.
+ void operator=(const Matrix<Form, Type>& mat);
- //- Assignment operator. Takes linear time.
- void operator=(const mType&);
+ //- Move assignment
+ void operator=(Matrix<Form, Type>&& mat);
//- Assignment to a block of another matrix
template<class MatrixType>
- void operator=(const ConstMatrixBlock<MatrixType>&);
+ void operator=(const ConstMatrixBlock<MatrixType>& Mb);
//- Assignment to a block of another matrix
template<class MatrixType>
- void operator=(const MatrixBlock<MatrixType>&);
+ void operator=(const MatrixBlock<MatrixType>& 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<Form, Type>& other);
+
+ //- Matrix subtraction
+ void operator-=(const Matrix<Form, Type>& other);
+
+ //- Matrix scalar multiplication
+ void operator*=(const scalar s);
+ //- Matrix scalar division
+ void operator/=(const scalar s);
- // IOstream operators
- //- Read Matrix from Istream, discarding contents of existing Matrix.
- friend Istream& operator>> <Form, Type>
- (
- Istream&,
- mType&
- );
+ // Random access iterator (non-const)
- //- Write Matrix to Ostream.
- friend Ostream& operator<< <Form, Type>
- (
- Ostream&,
- const mType&
- );
+ //- Return an iterator to begin traversing a Matrix
+ inline iterator begin();
+
+ //- 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<class Form, class Type>
-const Type& max(const Matrix<Form, Type>&);
+Istream& operator>>(Istream& is, Matrix<Form, Type>& mat);
+//- Write Matrix to Ostream, as per Matrix::writeMatrix() with
+//- default length, which is given by Detail::ListPolicy::short_length
template<class Form, class Type>
-const Type& min(const Matrix<Form, Type>&);
+Ostream& operator<<(Ostream& os, const Matrix<Form, Type>& mat);
+
+// Global Functions, Operators
+
+//- Find max value in the matrix
+template<class Form, class Type>
+const Type& max(const Matrix<Form, Type>& mat);
+
+//- Find min value in the matrix
+template<class Form, class Type>
+const Type& min(const Matrix<Form, Type>& mat);
+
+//- Find the min/max values of the matrix
+template<class Form, class Type>
+MinMax<Type> minMax(const Matrix<Form, Type>& mat);
+
+//- Matrix negation
template<class Form, class Type>
-Form operator-(const Matrix<Form, Type>&);
+Form operator-(const Matrix<Form, Type>& mat);
+//- Matrix addition
template<class Form, class Type>
Form operator+
(
- const Matrix<Form, Type>&,
- const Matrix<Form, Type>&
+ const Matrix<Form, Type>& A,
+ const Matrix<Form, Type>& B
);
+
+//- Matrix subtraction
template<class Form, class Type>
Form operator-
(
- const Matrix<Form, Type>&,
- const Matrix<Form, Type>&
+ const Matrix<Form, Type>& A,
+ const Matrix<Form, Type>& B
);
+
+//- Scalar multiplication of a matrix
template<class Form, class Type>
Form operator*
(
- const scalar,
- const Matrix<Form, Type>&
+ const scalar s,
+ const Matrix<Form, Type>& mat
);
+
+//- Scalar multiplication of a matrix
template<class Form, class Type>
Form operator*
(
- const Matrix<Form, Type>&,
- const scalar
+ const Matrix<Form, Type>& mat,
+ const scalar s
);
+
+//- Scalar division of a matrix
template<class Form, class Type>
Form operator/
(
- const Matrix<Form, Type>&,
- const scalar
+ const Matrix<Form, Type>& mat,
+ const scalar s
);
+
+//- Matrix-matrix multiplication
template<class Form1, class Form2, class Type>
typename typeOfInnerProduct<Type, Form1, Form2>::type
operator*
(
- const Matrix<Form1, Type>& a,
- const Matrix<Form2, Type>& b
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B
);
+
+//- Matrix-vector multiplication
template<class Form, class Type>
tmp<Field<Type>> operator*
(
- const Matrix<Form, Type>&,
- const Field<Type>&
+ const Matrix<Form, Type>& mat,
+ const Field<Type>& x
);
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index dfef8f70c5938c67c6f570c46304bbefdf2c5614..c264aed42ae5ff7a4d8df1f0d9e88e8c33b7369a 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -27,6 +27,20 @@ License
#include "MatrixBlock.H"
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+template<class Form, class Type>
+inline void Foam::Matrix<Form, Type>::doAlloc()
+{
+ const label len = size();
+
+ if (len)
+ {
+ v_ = new Type[len];
+ }
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Form, class Type>
@@ -39,8 +53,8 @@ inline Foam::Matrix<Form, Type>::Matrix()
template<class Form, class Type>
-inline Foam::autoPtr<Foam::Matrix<Form, Type>> Foam::Matrix<Form, Type>::
-clone() const
+inline Foam::autoPtr<Foam::Matrix<Form, Type>>
+Foam::Matrix<Form, Type>::clone() const
{
return autoPtr<Matrix<Form, Type>>::New(*this);
}
@@ -56,14 +70,14 @@ inline const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::m() const
+inline Foam::label Foam::Matrix<Form, Type>::m() const noexcept
{
return mRows_;
}
template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::n() const
+inline Foam::label Foam::Matrix<Form, Type>::n() const noexcept
{
return nCols_;
}
@@ -72,7 +86,14 @@ inline Foam::label Foam::Matrix<Form, Type>::n() const
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::size() const
{
- return mRows_*nCols_;
+ return mRows_ * nCols_;
+}
+
+
+template<class Form, class Type>
+inline bool Foam::Matrix<Form, Type>::empty() const noexcept
+{
+ return !mRows_ || !nCols_;
}
@@ -86,7 +107,7 @@ inline void Foam::Matrix<Form, Type>::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<Form, Type>::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<Form, Type>::checkj(const label j) const
template<class Form, class Type>
-inline const Type* Foam::Matrix<Form, Type>::v() const
+inline void Foam::Matrix<Form, Type>::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<class Form, class Type>
+inline bool Foam::Matrix<Form, Type>::uniform() const
+{
+ const label len = size();
+
+ if (len == 0)
+ {
+ return false;
+ }
+
+ for (label idx=1; idx<len; ++idx)
+ {
+ if (v_[0] != v_[idx])
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+
+template<class Form, class Type>
+inline const Type* Foam::Matrix<Form, Type>::cdata() const noexcept
{
return v_;
}
template<class Form, class Type>
-inline Type* Foam::Matrix<Form, Type>::v()
+inline Type* Foam::Matrix<Form, Type>::data() noexcept
{
return v_;
}
+template<class Form, class Type>
+inline const Type* Foam::Matrix<Form, Type>::rowData(const label irow) const
+{
+ #ifdef FULLDEBUG
+ checki(irow);
+ #endif
+ return (v_ + irow*nCols_);
+}
+
+
+template<class Form, class Type>
+inline Type* Foam::Matrix<Form, Type>::rowData(const label irow)
+{
+ #ifdef FULLDEBUG
+ checki(irow);
+ #endif
+ return (v_ + irow*nCols_);
+}
+
+
+template<class Form, class Type>
+inline const Type& Foam::Matrix<Form, Type>::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<class Form, class Type>
+inline Type& Foam::Matrix<Form, Type>::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<class Form, class Type>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
@@ -282,6 +388,13 @@ Foam::Matrix<Form, Type>::col
}
+template<class Form, class Type>
+inline void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
+{
+ resize(m, n);
+}
+
+
template<class Form, class Type>
void Foam::Matrix<Form, Type>::shallowResize(const label m, const label n)
{
@@ -290,47 +403,97 @@ void Foam::Matrix<Form, Type>::shallowResize(const label m, const label n)
}
+// * * * * * * * * * * * * * * * * * Iterators * * * * * * * * * * * * * * * //
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::iterator
+Foam::Matrix<Form, Type>::begin()
+{
+ return v_;
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::iterator
+Foam::Matrix<Form, Type>::end()
+{
+ return v_ + (mRows_ * nCols_);
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::const_iterator
+Foam::Matrix<Form, Type>::cbegin() const
+{
+ return v_;
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::const_iterator
+Foam::Matrix<Form, Type>::cend() const
+{
+ return v_ + (mRows_ * nCols_);
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::const_iterator
+Foam::Matrix<Form, Type>::begin() const
+{
+ return v_;
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::const_iterator
+Foam::Matrix<Form, Type>::end() const
+{
+ return v_ + (mRows_ * nCols_);
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Form, class Type>
inline const Type& Foam::Matrix<Form, Type>::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<class Form, class Type>
inline Type& Foam::Matrix<Form, Type>::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<class Form, class Type>
-inline const Type* Foam::Matrix<Form, Type>::operator[](const label i) const
+inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
{
- checki(i);
- return v_ + i*nCols_;
+ checki(irow);
+ return v_ + irow*nCols_;
}
template<class Form, class Type>
-inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
+inline Type* Foam::Matrix<Form, Type>::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 a39588225d2020d3603c5d616f999d6a82f8e1c8..3f9005c4025c93236a5725c2e396078780d4aaf9 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 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 <algorithm>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -45,26 +47,25 @@ Foam::Matrix<Form, Type>::Matrix(Istream& is)
template<class Form, class Type>
-Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
+bool Foam::Matrix<Form, Type>::readMatrix(Istream& is)
{
// Anull matrix
- M.clear();
+ clear();
is.fatalCheck(FUNCTION_NAME);
token firstToken(is);
- is.fatalCheck
- (
- "operator>>(Istream&, Matrix<Form, Type>&) : 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<Type>())
@@ -72,29 +73,22 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& 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<M.m(); i++)
+ // Loop over rows
+ for (label i=0; i < mRows_; ++i)
{
listDelimiter = is.readBeginList("MatrixRow");
- for (label j=0; j<M.n(); j++)
+ for (label j=0; j < nCols_; ++j)
{
- is >> v[k++];
+ is >> v_[idx++];
- is.fatalCheck
- (
- "operator>>(Istream&, Matrix<Form, Type>&) : "
- "reading entry"
- );
+ is.fatalCheck("readMatrix : reading reading entry");
}
is.readEndList("MatrixRow");
@@ -105,16 +99,9 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
Type element;
is >> element;
- is.fatalCheck
- (
- "operator>>(Istream&, Matrix<Form, Type>&) : "
- "reading the single entry"
- );
+ is.fatalCheck("readMatrix : reading the single entry");
- for (label i=0; i<mn; i++)
- {
- v[i] = element;
- }
+ std::fill(begin(), end(), element);
}
}
@@ -123,89 +110,72 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
}
else
{
- if (mn)
+ if (len)
{
- M.allocate();
- Type* v = M.v_;
+ is.read(reinterpret_cast<char*>(v_), len*sizeof(Type));
- is.read(reinterpret_cast<char*>(v), mn*sizeof(Type));
-
- is.fatalCheck
- (
- "operator>>(Istream&, Matrix<Form, Type>&) : "
- "reading the binary block"
- );
+ is.fatalCheck("readMatrix : reading the binary block");
}
}
- }
- else
- {
- FatalIOErrorInFunction(is)
- << "incorrect first token, expected <int>, found "
- << firstToken.info()
- << exit(FatalIOError);
+
+ return len;
}
- return is;
+
+ FatalIOErrorInFunction(is)
+ << "incorrect first token, expected <int>, found "
+ << firstToken.info()
+ << exit(FatalIOError);
+
+ return 0;
}
template<class Form, class Type>
-Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
+Foam::Ostream& Foam::Matrix<Form, Type>::writeMatrix
+(
+ Ostream& os,
+ const label shortLen
+) const
{
- label mn = M.mRows_*M.nCols_;
+ const Matrix<Form, Type>& 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<Type>())
{
- 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<Type>());
- if (uniform)
+ // Can the contents be considered 'uniform' (ie, identical)
+ if (len > 1 && contiguous<Type>() && mat.uniform())
{
- for (label i=0; i<mn; ++i)
- {
- if (v[i] != v[0])
- {
- uniform = false;
- break;
- }
- }
+ // Two or more entries, and all entries have identical values.
+ os << token::BEGIN_BLOCK << v[0] << token::END_BLOCK;
}
-
- if (uniform)
- {
- // Write start delimiter
- os << token::BEGIN_BLOCK;
-
- // Write contents
- os << v[0];
-
- // Write end delimiter
- os << token::END_BLOCK;
- }
- else if (mn < 10 && contiguous<Type>())
+ else if (len < shortLen && contiguous<Type>())
{
// Write start contents delimiter
os << token::BEGIN_LIST;
- label k = 0;
+ label idx = 0;
- // loop over rows
- for (label i=0; i<M.m(); i++)
+ // Loop over rows
+ for (label i=0; i < mat.m(); ++i)
{
os << token::BEGIN_LIST;
// Write row
- for (label j=0; j< M.n(); j++)
+ for (label j=0; j < mat.n(); ++j)
{
if (j) os << token::SPACE;
- os << v[k++];
+ os << v[idx++];
}
os << token::END_LIST;
@@ -219,17 +189,17 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& 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.m(); i++)
+ // Loop over rows
+ for (label i=0; i < mat.m(); ++i)
{
os << nl << token::BEGIN_LIST;
// Write row
- for (label j=0; j< M.n(); j++)
+ for (label j=0; j < mat.n(); ++j)
{
- os << nl << v[k++];
+ os << nl << v[idx++];
}
os << nl << token::END_LIST;
@@ -246,9 +216,16 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
}
else
{
- if (mn)
+ // Contents are binary and contiguous
+
+ if (len)
{
- os.write(reinterpret_cast<const char*>(M.v_), mn*sizeof(Type));
+ // write(...) includes surrounding start/end delimiters
+ os.write
+ (
+ reinterpret_cast<const char*>(mat.cdata()),
+ len*sizeof(Type)
+ );
}
}
@@ -257,4 +234,19 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
}
+template<class Form, class Type>
+Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& mat)
+{
+ mat.readMatrix(is);
+ return is;
+}
+
+
+template<class Form, class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& mat)
+{
+ return mat.writeMatrix(os, Detail::ListPolicy::short_length<Type>::value);
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index 75cd80a05271ad148c3a887d6d2d6d4fd50b5b68..1d754b5a198c328b272922ae09779825f0bd1fec 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<class MatrixType>
-Foam::QRMatrix<MatrixType>::QRMatrix(const MatrixType& M)
+Foam::QRMatrix<MatrixType>::QRMatrix(const MatrixType& mat)
{
- decompose(M);
+ decompose(mat);
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::decompose(const MatrixType& M)
+void Foam::QRMatrix<MatrixType>::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<MatrixType>::decompose(const MatrixType& M)
template<class MatrixType>
+template<template<typename> class ListContainer>
void Foam::QRMatrix<MatrixType>::solvex
(
- Field<cmptType>& x
+ ListContainer<cmptType>& x
) const
{
const label n = R_.n();
@@ -189,10 +190,11 @@ void Foam::QRMatrix<MatrixType>::solvex
template<class MatrixType>
void Foam::QRMatrix<MatrixType>::solve
(
- Field<cmptType>& x,
- const Field<cmptType>& source
+ List<cmptType>& x,
+ const UList<cmptType>& source
) const
{
+ // Assert (&x != &source) ?
const label m = Q_.m();
// x = Q_.T()*source;
@@ -213,15 +215,14 @@ template<class MatrixType>
Foam::tmp<Foam::Field<typename MatrixType::cmptType>>
Foam::QRMatrix<MatrixType>::solve
(
- const Field<cmptType>& source
+ const UList<cmptType>& source
) const
{
- tmp<Field<cmptType>> tx(new Field<cmptType>(Q_.m()));
- Field<cmptType>& x = tx.ref();
+ auto tresult(tmp<Field<cmptType>>::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 e925967646989300744743fb554185bd8bf5d0a8..b716cd195482b634f824c2285a749cf4411d9230 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 MatrixType>
class QRMatrix
{
-
public:
typedef typename MatrixType::cmptType cmptType;
typedef SquareMatrix<cmptType> 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<cmptType>& x) const;
+ //- the appropriate transformed source (e.g. Q.T()*source)
+ //- and return the solution in x
+ template<template<typename> class ListContainer>
+ void solvex(ListContainer<cmptType>& 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<cmptType>& x, const Field<cmptType>& source) const;
+ //- and return the solution in the Field argument x
+ void solve(List<cmptType>& x, const UList<cmptType>& source) const;
//- Solve the linear system with the given source
- // returning the solution
- tmp<Field<cmptType>> solve(const Field<cmptType>& source) const;
+ //- and return the solution
+ tmp<Field<cmptType>> solve(const UList<cmptType>& 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 256c22bac40abd8f1255a67d90ad749b0b0f235d..c1d6a0826bb602d461a7718127ac093ecead2167 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 924615f9eb905eebdf9f2dfa259ecbd94525d201..ae49b677e51d57c4eec8edb09493b923e183a8a2 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
@@ -77,27 +77,30 @@ public:
inline RectangularMatrix(const MatrixBlock<MatrixType>&);
//- 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<Type>&);
+ inline RectangularMatrix(const SquareMatrix<Type>& mat);
//- Construct from Istream.
- inline RectangularMatrix(Istream&);
+ inline RectangularMatrix(Istream& is);
//- Clone
inline autoPtr<RectangularMatrix<Type>> 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 6d40febe009df0450fd5c1f06f0c9d7a53552963..f45b3b797803ee4625265b241b29c0dcd74dd105 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
@@ -84,20 +84,20 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const label m,
const label n,
- const Type& t
+ const Type& val
)
:
- Matrix<RectangularMatrix<Type>, Type>(m, n, t)
+ Matrix<RectangularMatrix<Type>, Type>(m, n, val)
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
- const SquareMatrix<Type>& SM
+ const SquareMatrix<Type>& mat
)
:
- Matrix<RectangularMatrix<Type>, Type>(SM)
+ Matrix<RectangularMatrix<Type>, Type>(mat)
{}
@@ -125,6 +125,13 @@ void Foam::RectangularMatrix<Type>::operator=(const zero)
}
+template<class Type>
+void Foam::RectangularMatrix<Type>::operator=(const Type& val)
+{
+ Matrix<RectangularMatrix<Type>, Type>::operator=(val);
+}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index 576f95db793d7cf331c23d6a48a5dea1ab20d7be..f1819a9c9cb46385ddacb248f52b64ef4589a733 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 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 Type>
-class RectangularMatrix;
+// Forward Declarations
+template<class Type> class RectangularMatrix;
/*---------------------------------------------------------------------------*\
@@ -82,23 +80,23 @@ public:
inline SquareMatrix(const MatrixBlock<MatrixType>&);
//- 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<Type>);
//- 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<Type>&);
//- 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<Type>);
};
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index c65f617edb309896f1a0720edef7911170dce38f..97c61c892372a0da020ec6d1010dabb497572ec7 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 | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -98,42 +98,42 @@ template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label n,
- const Identity<Type>
+ const Type& val
)
:
- Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
-{
- for (label i=0; i<n; i++)
- {
- this->operator()(i, i) = Type(I);
- }
-}
+ Matrix<SquareMatrix<Type>, Type>(n, n, val)
+{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label n,
- const Type& t
+ const Identity<Type>
)
:
- Matrix<SquareMatrix<Type>, Type>(n, n, t)
-{}
+ Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
+{
+ for (label i=0; i<n; i++)
+ {
+ this->operator()(i, i) = Type(I);
+ }
+}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
- const RectangularMatrix<Type>& RM
+ const RectangularMatrix<Type>& mat
)
:
- Matrix<SquareMatrix<Type>, Type>(RM)
+ Matrix<SquareMatrix<Type>, 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<Type>::clone() const
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+template<class Type>
+inline void Foam::SquareMatrix<Type>::resize(const label m)
+{
+ Matrix<SquareMatrix<Type>, Type>::resize(m, m);
+}
+
+
template<class Type>
inline void Foam::SquareMatrix<Type>::setSize(const label m)
{
- Matrix<SquareMatrix<Type>, Type>::setSize(m, m);
+ Matrix<SquareMatrix<Type>, Type>::resize(m, m);
}
@@ -179,10 +186,18 @@ void Foam::SquareMatrix<Type>::operator=(const zero)
}
+template<class Type>
+void Foam::SquareMatrix<Type>::operator=(const Type& val)
+{
+ Matrix<SquareMatrix<Type>, Type>::operator=(val);
+}
+
+
template<class Type>
void Foam::SquareMatrix<Type>::operator=(const Identity<Type>)
{
Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
+
for (label i=0; i<this->n(); i++)
{
this->operator()(i, i) = Type(I);
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
index 267da6b327f5114a270240e7401f85d756b0075f..f4bf7c9c41aaf24574b79b8528877e8f464d9bd7 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<Type>);
- //- 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<class Type>
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
index 7dcc04e2b7527b60aee5cb5f511c1a160b96cdf2..a2f8dc01092d2bc367084a4f1f90c581350faedf 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
@@ -56,27 +56,27 @@ template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label n,
- const Identity<Type>
+ const Type& val
)
:
- Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
-{
- for (label i=0; i<n; i++)
- {
- this->operator()(i, i) = pTraits<Type>::one;
- }
-}
+ Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, val)
+{}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label n,
- const Type& t
+ const Identity<Type>
)
:
- Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, t)
-{}
+ Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
+{
+ for (label i=0; i<n; i++)
+ {
+ this->operator()(i, i) = pTraits<Type>::one;
+ }
+}
template<class Type>