Commit 1a60bb4c authored by Mark Olesen's avatar Mark Olesen

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
parent a9688d07
......@@ -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));
......
......@@ -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;
......
......@@ -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);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -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;
}
......
......@@ -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;
};
......
......@@ -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;
......
......@@ -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;
}
......
This diff is collapsed.
This diff is collapsed.
......@@ -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];
}
}