zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: improvements, modernization of matrix containers (#1220)

Browse Source 
- 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
This commit is contained in:
Mark Olesen
2019-05-22 12:18:31 +01:00
committed by Andrew Heather
parent f8a70115fd
commit 061eb53fb5
20 changed files with 1121 additions and 730 deletions
Download Patch File Download Diff File Expand all files Collapse all files

114
applications/test/Matrix/Test-Matrix.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -40,44 +40,82 @@ using namespace Foam;
int main(int argc, char *argv[])
{
SquareMatrix<scalar> hmm(3);
hmm(0, 0) = -3.0;
hmm(0, 1) = 10.0;
hmm(0, 2) = -4.0;
hmm(1, 0) = 2.0;
hmm(1, 1) = 3.0;
hmm(1, 2) = 10.0;
hmm(2, 0) = 2.0;
hmm(2, 1) = 6.0;
hmm(2, 2) = 1.0;
//Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
Info<< max(hmm) << endl;
Info<< min(hmm) << endl;
SquareMatrix<scalar> hmm2(3, I);
hmm = hmm2;
Info<< hmm << endl;
//SquareMatrix<scalar> hmm3(Sin);
//Info<< hmm3 << endl;
SquareMatrix<scalar> hmm4;
hmm4 = hmm2;
Info<< hmm4 << endl;
SquareMatrix<scalar> hmm5;
hmm4 = hmm5;
Info<< hmm5 << endl;
// 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;
Info<< "min/max: " << minMax(square1) << nl;
square1 = 100;
Info<< "matrix: " << square1
<< " begin: " << long(square1.cdata()) << nl;
SquareMatrix<scalar> square2(3, I);
square1 = square2;
Info<< nl << "After copy assign from Identity:" << nl << square1 << nl;
square1 += 1.25*square2;
Info<< nl << "After +=" << nl << square1 << nl;
square1 -= square2.T();
Info<< nl << "After -=" << nl << square1 << nl;
square1 *= 10;
Info<< nl << "After *=" << nl << square1 << nl;
square1 /= 8;
Info<< nl << "After /=" << nl << square1 << nl;
SquareMatrix<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));

74
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
View File

@ -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
@ -36,46 +36,55 @@ inline Foam::DiagonalMatrix<Type>::DiagonalMatrix()
{}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n)
:
List<Type>(n)
{}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n, const zero)
:
List<Type>(n, Zero)
{}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n, const Type& val)
:
List<Type>(n, val)
{}
template<class Type>
template<class Form>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& mat)
:
List<Type>(min(a.m(), a.n()))
List<Type>(min(mat.m(), mat.n()))
{
forAll(*this, i)
label i = 0;
for (Type& val : *this)
{
this->operator[](i) = a(i, i);
val = mat(i, i);
++i;
}
}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size)
:
List<Type>(size)
{}
template<class Type>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size, const Type& val)
:
List<Type>(size, val)
{}
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;

31
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -27,8 +27,7 @@ Class
Foam::DiagonalMatrix
Description
DiagonalMatrix<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,8 +44,7 @@ SourceFiles
namespace Foam
{
// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * //
// Forward Declarations
template<class Form, class Type> class Matrix;
/*---------------------------------------------------------------------------*\
@ -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);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

40
src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
View File

@ -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;
}

19
src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H
View File

@ -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;
};

4
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
View File

@ -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;

35
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -26,15 +26,15 @@ License
\*---------------------------------------------------------------------------*/
#include "LUscalarMatrix.H"
#include "SubField.H"
#include "SubList.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<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
(
@ -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;
}

575
src/OpenFOAM/matrices/Matrix/Matrix.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -26,18 +26,8 @@ License
\*---------------------------------------------------------------------------*/
#include "Matrix.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<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();
doAlloc();
std::fill(begin(), end(), val);
}
allocate();
if (v_)
template<class Form, class Type>
Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& mat)
:
mRows_(mat.mRows_),
nCols_(mat.nCols_),
v_(nullptr)
{
const label mn = size();
for (label i=0; i<mn; i++)
if (mat.cdata())
{
v_[i] = s;
}
doAlloc();
std::copy(mat.cbegin(), mat.cend(), v_);
}
}
template<class Form, class Type>
Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& M)
Foam::Matrix<Form, Type>::Matrix(Matrix<Form, Type>&& mat)
:
mRows_(M.mRows_),
nCols_(M.nCols_),
v_(nullptr)
mRows_(mat.mRows_),
nCols_(mat.nCols_),
v_(mat.v_)
{
if (M.v_)
{
allocate();
const label mn = size();
for (label i=0; i<mn; i++)
{
v_[i] = M.v_[i];
}
}
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,11 +130,11 @@ 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);
}
@ -186,11 +152,11 @@ 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);
}
@ -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;
nCols_ = M.nCols_;
M.nCols_ = 0;
v_ = M.v_;
M.v_ = nullptr;
return list;
}
template<class Form, class Type>
void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
void Foam::Matrix<Form, Type>::swap(Matrix<Form, Type>& mat)
{
mType newMatrix(m, n, Zero);
Foam::Swap(mRows_, mat.mRows_);
Foam::Swap(nCols_, mat.nCols_);
Foam::Swap(v_, mat.v_);
}
label minM = min(m, mRows_);
label minN = min(n, nCols_);
for (label i=0; i<minM; i++)
template<class Form, class Type>
void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& mat)
{
for (label j=0; j<minN; j++)
clear();
mRows_ = mat.mRows_;
nCols_ = mat.nCols_;
v_ = mat.v_;
mat.mRows_ = 0;
mat.nCols_ = 0;
mat.v_ = nullptr;
}
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 < 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++)
std::copy(mat.cbegin(), mat.cend(), v_);
}
}
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator=(Matrix<Form, Type>&& mat)
{
v_[i] = M.v_[i];
}
if (this == &mat)
{
FatalErrorInFunction
<< "Attempted assignment to self"
<< abort(FatalError);
}
this->transfer(mat);
}
@ -318,9 +326,9 @@ 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);
}
@ -335,9 +343,9 @@ 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);
}
@ -346,29 +354,97 @@ void Foam::Matrix<Form, Type>::operator=
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);
}
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& other)
{
if (this == &other)
{
FatalErrorInFunction
<< "Attempted addition to self"
<< abort(FatalError);
}
if (m() != other.m() || n() != other.n())
{
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];
}
}
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);
}
if (m() != other.m() || n() != other.n())
{
FatalErrorInFunction
<< "Attempt to subtract matrices with different sizes: ("
<< m() << ", " << n() << ") != ("
<< other.m() << ", " << other.n() << ')' << nl
<< abort(FatalError);
}
Type* out = this->data();
const Type* in = other.cdata();
const label len = this->size();
for (label idx=0; idx < len; ++idx)
{
out[idx] -= in[idx];
}
}
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator*=(const scalar s)
{
for (Type& val : *this)
{
val *= s;
}
}
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator/=(const scalar s)
{
for (Type& val : *this)
{
val /= s;
}
}
@ -376,89 +452,61 @@ void Foam::Matrix<Form, Type>::operator=(const zero)
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
template<class Form, class Type>
const Type& Foam::max(const Matrix<Form, Type>& M)
const Type& Foam::max(const Matrix<Form, Type>& mat)
{
const label mn = M.size();
if (mn)
{
label curMaxI = 0;
const Type* Mv = M.v();
for (label i=1; i<mn; i++)
{
if (Mv[i] > Mv[curMaxI])
{
curMaxI = i;
}
}
return Mv[curMaxI];
}
else
if (mat.empty())
{
FatalErrorInFunction
<< "Matrix is empty"
<< abort(FatalError);
// Return in error to keep compiler happy
return M(0, 0);
<< "Matrix is empty" << abort(FatalError);
}
return *(std::max_element(mat.cbegin(), mat.cend()));
}
template<class Form, class Type>
const Type& Foam::min(const Matrix<Form, Type>& M)
const Type& Foam::min(const Matrix<Form, Type>& mat)
{
const label mn = M.size();
if (mn)
{
label curMinI = 0;
const Type* Mv = M.v();
for (label i=1; i<mn; i++)
{
if (Mv[i] < Mv[curMinI])
{
curMinI = i;
}
}
return Mv[curMinI];
}
else
if (mat.empty())
{
FatalErrorInFunction
<< "Matrix is empty"
<< abort(FatalError);
// Return in error to keep compiler happy
return M(0, 0);
<< "Matrix is empty" << abort(FatalError);
}
return *(std::min_element(mat.cbegin(), mat.cend()));
}
template<class Form, class Type>
Foam::MinMax<Type> Foam::minMax(const Matrix<Form, Type>& mat)
{
MinMax<Type> result;
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());
Form result(mat.m(), mat.n());
if (M.m() && M.n())
{
Type* nMv = nM.v();
const Type* Mv = M.v();
std::transform
(
mat.cbegin(),
mat.cend(),
result.begin(),
std::negate<Type>()
);
const label mn = M.size();
for (label i=0; i<mn; i++)
{
nMv[i] = -Mv[i];
}
}
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 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;
@ -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());
if (M.m() && M.n())
{
Type* sMv = sM.v();
const Type* Mv = M.v();
const label len = mat.size();
const label mn = M.size();
for (label i=0; i<mn; i++)
if (len)
{
sMv[i] = s*Mv[i];
Type* out = result.data();
const Type* in = mat.cdata();
for (label idx=0; idx < len; ++idx)
{
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())
{
Type* sMv = sM.v();
const Type* Mv = M.v();
const label len = mat.size();
const label mn = M.size();
for (label i=0; i<mn; i++)
if (len)
{
sMv[i] = Mv[i]*s;
Type* out = result.data();
const Type* in = mat.cdata();
for (label idx=0; idx < len; ++idx)
{
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());
if (M.m() && M.n())
{
Type* sMv = sM.v();
const Type* Mv = M.v();
const label len = mat.size();
const label mn = M.size();
for (label i=0; i<mn; i++)
if (len)
{
sMv[i] = Mv[i]/s;
Type* out = result.data();
const Type* in = mat.cdata();
for (label idx=0; idx < len; ++idx)
{
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,11 +659,11 @@ 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);
}
@ -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;
}

321
src/OpenFOAM/matrices/Matrix/Matrix.H
View File

@ -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
@ -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 constructor.
Matrix(const mType&);
//- Copy construct
Matrix(const Matrix<Form, Type>& mat);
//- 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;
@ -154,23 +156,45 @@ public:
// 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 element vector of the constant Matrix
inline const Type* v() const;
//- Return true if the matrix is empty (ie, size() is zero)
inline bool empty() const noexcept;
//- Return element vector of the Matrix
inline Type* v();
//- Return const pointer to the first data element, which can also
//- be used to address into the matrix contents
inline const Type* cdata() const noexcept;
//- Return pointer to the first data element, which can also
//- be used to address into the matrix contents
inline Type* data() noexcept;
//- Return const pointer to data in the specified row.
// Subscript checking only with FULLDEBUG
inline const Type* rowData(const label irow) const;
//- Return pointer to data in the specified row.
// Subscript checking only with FULLDEBUG
inline Type* rowData(const label irow);
//- Linear addressing const element access
// Subscript checking only with FULLDEBUG
inline const Type& at(const label idx) const;
//- Linear addressing element access
// Subscript checking only with FULLDEBUG
inline Type& at(const label idx);
// Block access
// Block Access
inline ConstMatrixBlock<mType> block
(
@ -229,14 +253,19 @@ public:
const label nStart
);
// Check
//- Check index i is within valid range (0 ... m-1).
inline void checki(const label i) const;
inline void checki(const label irow) const;
//- Check index j is within valid range (0 ... n-1).
inline void checkj(const label j) const;
inline void checkj(const label jcol) const;
//- Check that dimensions are positive, non-zero
inline void checkSize() const;
//- True if all entries have identical values, and matrix is non-empty
inline bool uniform() const;
// Edit
@ -244,12 +273,22 @@ public:
//- 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);
//- 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);
@ -259,115 +298,201 @@ public:
Form T() const;
// Member operators
// Member Operators
//- Return subscript-checked row of Matrix.
inline Type* operator[](const label);
//- Return const pointer to data in the specified row - rowData().
// Subscript checking only with FULLDEBUG
inline const Type* operator[](const label irow) const;
//- Return subscript-checked row of constant Matrix.
inline const Type* operator[](const label) const;
//- Return pointer to data in the specified row - rowData().
// Subscript checking only with FULLDEBUG
inline Type* operator[](const label irow);
//- (i, j) const element access operator
inline const Type& operator()(const label i, const label j) const;
// Subscript checking only with FULLDEBUG
inline const Type& operator()(const label irow, const label jcol) const;
//- (i, j) element access operator
inline Type& operator()(const label i, const label j);
// Subscript checking only with FULLDEBUG
inline Type& operator()(const label irow, const label jcol);
//- Assignment operator. Takes linear time.
void operator=(const mType&);
//- Copy assignment. Takes linear time.
void operator=(const Matrix<Form, Type>& mat);
//- 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
// Random access iterator (non-const)
//- Read Matrix from Istream, discarding contents of existing Matrix.
friend Istream& operator>> <Form, Type>
(
Istream&,
mType&
);
//- Return an iterator to begin traversing a Matrix
inline iterator begin();
//- Write Matrix to Ostream.
friend Ostream& operator<< <Form, Type>
(
Ostream&,
const mType&
);
//- Return an iterator to end traversing a Matrix
inline iterator end();
// Random access iterator (const)
//- Return const_iterator to begin traversing a constant Matrix
inline const_iterator cbegin() const;
//- Return const_iterator to end traversing a constant Matrix
inline const_iterator cend() const;
//- Return const_iterator to begin traversing a constant Matrix
inline const_iterator begin() const;
//- Return const_iterator to end traversing a constant Matrix
inline const_iterator end() const;
// IO
//- Read Matrix from Istream, discarding existing contents.
bool readMatrix(Istream& is);
//- Write Matrix, with line-breaks in ASCII when length exceeds
//- shortLen.
// Using '0' suppresses line-breaks entirely.
Ostream& writeMatrix(Ostream& os, const label shortLen=0) const;
// Housekeeping
//- Deprecated(2019-04) raw data pointer, const access
// \deprecated(2019-04) - use cdata() method
const Type* FOAM_DEPRECATED_FOR(2019-04, "cdata() method") v() const
{
return v_;
}
//- Deprecated(2019-04) raw data pointer, non-const access
// \deprecated(2019-04) - use data() method
Type* FOAM_DEPRECATED_FOR(2019-04, "data() method") v()
{
return v_;
}
};
// Global functions and operators
// IOstream Operators
//- Read Matrix from Istream, discarding contents of existing Matrix.
template<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>
Form operator-(const Matrix<Form, 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>& 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
);

213
src/OpenFOAM/matrices/Matrix/MatrixI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -27,6 +27,20 @@ License
#include "MatrixBlock.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<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_;
}
@ -76,6 +90,13 @@ inline Foam::label Foam::Matrix<Form, Type>::size() const
}
template<class Form, class Type>
inline bool Foam::Matrix<Form, Type>::empty() const noexcept
{
return !mRows_ || !nCols_;
}
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checki(const label i) const
{
@ -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_;
}

170
src/OpenFOAM/matrices/Matrix/MatrixIO.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -30,6 +30,8 @@ License
#include "Ostream.H"
#include "token.H"
#include "contiguous.H"
#include "ListPolicy.H"
#include <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("readMatrix : reading the binary block");
}
}
is.fatalCheck
(
"operator>>(Istream&, Matrix<Form, Type>&) : "
"reading the binary block"
);
return len;
}
}
}
else
{
FatalIOErrorInFunction(is)
<< "incorrect first token, expected <int>, found "
<< firstToken.info()
<< exit(FatalIOError);
}
return is;
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);
}
// ************************************************************************* //

31
src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
View File

@ -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;
}

25
src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
View File

@ -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;
};

2
src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
View File

@ -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

15
src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
View File

@ -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);
};

15
src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
View File

@ -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

28
src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -47,10 +47,8 @@ SourceFiles
namespace Foam
{
// Forward declaration of friend functions and operators
template<class 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.
@ -112,6 +110,9 @@ public:
// Edit
//- Resize the matrix preserving the elements
inline void resize(const label m);
//- Resize the matrix preserving the elements
inline void setSize(const label m);
@ -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>);
};

49
src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -94,6 +94,17 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label n,
const Type& val
)
:
Matrix<SquareMatrix<Type>, Type>(n, n, val)
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
@ -113,27 +124,16 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label n,
const Type& t
const RectangularMatrix<Type>& mat
)
:
Matrix<SquareMatrix<Type>, Type>(n, n, t)
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const RectangularMatrix<Type>& RM
)
:
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);

12
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
View File

@ -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>

24
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -52,6 +52,17 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
{}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label n,
const Type& val
)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, val)
{}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
@ -68,17 +79,6 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label n,
const Type& t
)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, t)
{}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(Istream& is)
: