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: harmonize matrix constructors (#1220)

Browse Source 
- generalize identity matrix constructors for non-scalar types

- add constructors using labelPair for the row/column sizing information.
  For a SquareMatrix, this provides an unambiguous parameter resolution.

- reuse assignment operators

STYLE: adjust matrix comments
This commit is contained in:
Mark Olesen
2019-05-29 09:50:46 +02:00
committed by Andrew Heather
parent 2bdcd5b80d
commit 96d0a8f2af
18 changed files with 522 additions and 185 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -33,8 +33,63 @@ License
#include "tensor.H"
#include "IFstream.H"
#include <algorithm>
using namespace Foam;
// Copy values into matrix
template<class Form, class Type>
void assignMatrix
(
Matrix<Form, Type>& mat,
std::initializer_list<typename Matrix<Form, Type>::cmptType> list
)
{
const label nargs = list.size();
if (nargs != mat.size())
{
FatalErrorInFunction
<< "Mismatch in matrix dimension ("
<< mat.m() << ", "
<< mat.n() << ") and number of args (" << nargs << ')' << nl
<< exit(FatalError);
}
std::copy(list.begin(), list.end(), mat.begin());
}
// Create matrix with values
template<class MatrixType>
MatrixType makeMatrix
(
const labelPair& dims,
std::initializer_list<typename MatrixType::cmptType> list
)
{
MatrixType mat(dims);
assignMatrix(mat, list);
return mat;
}
// Create matrix with values
template<class MatrixType, Foam::label nRows, Foam::label nCols>
MatrixType makeMatrix
(
std::initializer_list<typename MatrixType::cmptType> list
)
{
MatrixType mat(labelPair(nRows, nCols));
assignMatrix(mat, list);
return mat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
@ -46,18 +101,48 @@ int main(int argc, char *argv[])
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;
assignMatrix
(
square1,
{
-3.0, 10.0, -4.0,
2.0, 3.0, 10.0,
2.0, 6.0, 1.0
}
);
Info<< "matrix: " << square1
<< " begin: " << long(square1.cdata()) << nl;
<< " begin: " << uintptr_t(square1.cdata()) << nl;
// Test makeMatrix
{
auto square1b
(
makeMatrix<scalarSquareMatrix>
(
{3, 3},
{
-3.0, 10.0, -4.0,
2.0, 3.0, 10.0,
2.0, 6.0, 1.0
}
)
);
Info<< "makeMatrix: " << square1b << nl;
auto square1c
(
makeMatrix<scalarSquareMatrix, 3, 3>
({
-3.0, 10.0, -4.0,
2.0, 3.0, 10.0,
2.0, 6.0, 1.0
})
);
Info<< "makeMatrix: " << square1c << nl;
}
//Info<< square1 - 2.0*(-square1) << nl;
Info<< "min:" << min(square1) << " max:" << max(square1) << nl;
@ -68,7 +153,7 @@ int main(int argc, char *argv[])
List<scalar> stole(square1.release());
Info<< "matrix: " << square1
<< " begin: " << long(square1.cdata()) << nl;
<< " begin: " << uintptr_t(square1.cdata()) << nl;
Info<< "List: " << stole << nl;
@ -77,7 +162,7 @@ int main(int argc, char *argv[])
square1 = 100;
Info<< "matrix: " << square1
<< " begin: " << long(square1.cdata()) << nl;
<< " begin: " << uintptr_t(square1.cdata()) << nl;
SquareMatrix<scalar> square2(3, I);
@ -116,10 +201,13 @@ int main(int argc, char *argv[])
{
Info<< nl << "Test RectangularMatrix" << nl;
RectangularMatrix<scalar> rm1(5, 6, 3.1);
RectangularMatrix<scalar> rm1({5, 6}, 3.1);
rm1(0, 1) = 4.5;
RectangularMatrix<scalar> rm1b(rm1.block(2, 2, 0, 0));
Info<< "rm1b = " << rm1b << endl;
Info // << "Full matrix " << rm1 << nl
<< "block = " << rm1b << endl;
}
{
@ -143,7 +231,7 @@ int main(int argc, char *argv[])
Info<< "Inverse = " << invDecomposed(symmMatrix2) << endl;
Info<< "Determinant = " << detDecomposed(symmMatrix2) << endl;
scalarDiagonalMatrix rhs(3, 0);
scalarDiagonalMatrix rhs(3, Zero);
rhs[0] = 1;
rhs[1] = 2;
rhs[2] = 3;

2
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
View File

@ -145,7 +145,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
else
{
label nCells = ldum.lduAddr().size();
scalarSquareMatrix m(nCells, 0.0);
scalarSquareMatrix m(nCells, Zero);
transfer(m);
convert(ldum, interfaceCoeffs, interfaces);
}

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

@ -107,7 +107,7 @@ public:
//- Construct from lduMatrix and perform LU decomposition
LUscalarMatrix
(
const lduMatrix&,
const lduMatrix& ldum,
const FieldField<Field, scalar>& interfaceCoeffs,
const lduInterfaceFieldPtrsList& interfaces
);

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

@ -586,7 +586,7 @@ Form Foam::operator-(const Matrix<Form, Type>& mat)
template<class Form, class Type>
Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
{
if (A.m() != B.m())
if (A.m() != B.m() || A.n() != B.n())
{
FatalErrorInFunction
<< "Attempt to add matrices with different sizes: ("
@ -615,7 +615,7 @@ Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
template<class Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
{
if (A.m() != B.m())
if (A.m() != B.m() || A.n() != B.n())
{
FatalErrorInFunction
<< "Attempt to subtract matrices with different sizes: ("

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

@ -50,6 +50,7 @@ SourceFiles
#define Matrix_H
#include "uLabel.H"
#include "Pair.H"
#include "Field.H"
#include "autoPtr.H"
@ -120,20 +121,31 @@ public:
// Constructors
//- Null constructor.
//- Construct null
inline Matrix();
//- Construct given number of rows and columns.
//- Construct given number of rows/columns
Matrix(const label m, const label n);
//- Construct with given number of rows and columns
//- Construct with given number of rows/columns
//- initializing all elements to zero
Matrix(const label m, const label n, const zero);
//- Construct with given number of rows and columns
//- Construct with given number of rows/columns
//- initializing all elements to the given value
Matrix(const label m, const label n, const Type& val);
//- Construct given number of rows/columns
inline explicit Matrix(const labelPair& dims);
//- Construct given number of rows/columns.
//- initializing all elements to zero
inline Matrix(const labelPair& dims, const zero);
//- Construct with given number of rows/columns
//- initializing all elements to the given value
inline Matrix(const labelPair& dims, const Type& val);
//- Copy construct
Matrix(const Matrix<Form, Type>& mat);
@ -176,6 +188,9 @@ public:
//- Return the number of elements in matrix (m*n)
inline label size() const;
//- Return row/column sizes
inline labelPair sizes() const;
//- Return true if the matrix is empty (ie, size() is zero)
inline bool empty() const noexcept;

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

@ -52,6 +52,27 @@ inline Foam::Matrix<Form, Type>::Matrix()
{}
template<class Form, class Type>
inline Foam::Matrix<Form, Type>::Matrix(const labelPair& dims)
:
Matrix<Form, Type>(dims.first(), dims.second())
{}
template<class Form, class Type>
inline Foam::Matrix<Form, Type>::Matrix(const labelPair& dims, const zero)
:
Matrix<Form, Type>(dims.first(), dims.second(), Zero)
{}
template<class Form, class Type>
inline Foam::Matrix<Form, Type>::Matrix(const labelPair& dims, const Type& val)
:
Matrix<Form, Type>(dims.first(), dims.second(), val)
{}
template<class Form, class Type>
inline Foam::autoPtr<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::clone() const
@ -90,6 +111,13 @@ inline Foam::label Foam::Matrix<Form, Type>::size() const
}
template<class Form, class Type>
inline Foam::labelPair Foam::Matrix<Form, Type>::sizes() const
{
return labelPair(mRows_, nCols_);
}
template<class Form, class Type>
inline bool Foam::Matrix<Form, Type>::empty() const noexcept
{

12
src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.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
@ -60,7 +60,7 @@ namespace Foam
template<class MatrixType>
class ConstMatrixBlock
{
// Private data
// Private Data
//- Const reference to the parent matrix
const MatrixType& matrix_;
@ -98,6 +98,9 @@ public:
//- Return the number of columns in the block
inline label n() const;
//- Return row/column sizes
inline labelPair sizes() const;
//- (i, j) const element access operator
inline const cmptType& operator()
(
@ -117,7 +120,7 @@ public:
template<class MatrixType>
class MatrixBlock
{
// Private data
// Private Data
//- Reference to the parent matrix
MatrixType& matrix_;
@ -155,6 +158,9 @@ public:
//- Return the number of columns in the block
inline label n() const;
//- Return row/column sizes
inline labelPair sizes() const;
//- (i, j) const element access operator
inline const cmptType& operator()
(

16
src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.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
@ -119,6 +119,20 @@ inline Foam::label Foam::MatrixBlock<MatrixType>::n() const
}
template<class MatrixType>
inline Foam::labelPair Foam::ConstMatrixBlock<MatrixType>::sizes() const
{
return labelPair(mRows_, nCols_);
}
template<class MatrixType>
inline Foam::labelPair Foam::MatrixBlock<MatrixType>::sizes() const
{
return labelPair(mRows_, nCols_);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class MatrixType>

54
src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.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
@ -49,7 +49,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
Class RectangularMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
@ -62,33 +62,47 @@ public:
// Constructors
//- Null constructor.
//- Construct null
inline RectangularMatrix();
//- Construct given number of rows and columns,
//- Construct a square matrix (rows == columns)
inline explicit RectangularMatrix(const label n);
//- Construct given number of rows/columns
inline RectangularMatrix(const label m, const label n);
//- Construct from a block of another matrix
template<class MatrixType>
inline RectangularMatrix(const ConstMatrixBlock<MatrixType>&);
//- Construct from a block of another matrix
template<class MatrixType>
inline RectangularMatrix(const MatrixBlock<MatrixType>&);
//- Construct with given number of rows and columns
//- Construct given number of rows/columns
//- 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.
//- Construct given number of rows/columns
//- initializing all elements to the given value
inline RectangularMatrix(const label m, const label n, const Type& val);
//- Construct given number of rows/columns
inline explicit RectangularMatrix(const labelPair& dims);
//- Construct given number of rows/columns
//- initializing all elements to zero
inline RectangularMatrix(const labelPair& dims, const zero);
//- Construct given number of rows/columns
//- initializing all elements to the given value
inline RectangularMatrix(const labelPair& dims, const Type& val);
//- Construct from a block of another matrix
template<class MatrixType>
inline RectangularMatrix(const ConstMatrixBlock<MatrixType>& mat);
//- Construct from a block of another matrix
template<class MatrixType>
inline RectangularMatrix(const MatrixBlock<MatrixType>& mat);
//- Construct as copy of a square matrix
inline RectangularMatrix(const SquareMatrix<Type>& mat);
//- Construct from Istream.
inline RectangularMatrix(Istream& is);
inline explicit RectangularMatrix(Istream& is);
//- Clone
inline autoPtr<RectangularMatrix<Type>> clone() const;
@ -96,11 +110,11 @@ public:
// Member Operators
//- Assignment of all elements to zero
void operator=(const zero);
//- Assign all elements to zero
inline void operator=(const zero);
//- Assignment of all elements to the given value
void operator=(const Type& val);
//- Assign all elements to value
inline void operator=(const Type& val);
};

77
src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.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
@ -46,24 +46,12 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
template<class Type>
template<class MatrixType>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const ConstMatrixBlock<MatrixType>& block
const label n
)
:
Matrix<RectangularMatrix<Type>, Type>(block)
{}
template<class Type>
template<class MatrixType>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const MatrixBlock<MatrixType>& block
)
:
Matrix<RectangularMatrix<Type>, Type>(block)
Matrix<RectangularMatrix<Type>, Type>(n, n)
{}
@ -91,6 +79,60 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const labelPair& dims
)
:
RectangularMatrix<Type>(dims.first(), dims.second())
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const labelPair& dims,
const zero
)
:
RectangularMatrix<Type>(dims.first(), dims.second(), Zero)
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const labelPair& dims,
const Type& val
)
:
RectangularMatrix<Type>(dims.first(), dims.second(), val)
{}
template<class Type>
template<class MatrixType>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const ConstMatrixBlock<MatrixType>& mat
)
:
Matrix<RectangularMatrix<Type>, Type>(mat)
{}
template<class Type>
template<class MatrixType>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const MatrixBlock<MatrixType>& mat
)
:
Matrix<RectangularMatrix<Type>, Type>(mat)
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
@ -119,14 +161,14 @@ Foam::RectangularMatrix<Type>::clone() const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
void Foam::RectangularMatrix<Type>::operator=(const zero)
inline void Foam::RectangularMatrix<Type>::operator=(const zero)
{
Matrix<RectangularMatrix<Type>, Type>::operator=(Zero);
}
template<class Type>
void Foam::RectangularMatrix<Type>::operator=(const Type& val)
inline void Foam::RectangularMatrix<Type>::operator=(const Type& val)
{
Matrix<RectangularMatrix<Type>, Type>::operator=(val);
}
@ -164,4 +206,5 @@ inline Foam::RectangularMatrix<Type> outer
} // End namespace Foam
// ************************************************************************* //

15
src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
View File

@ -72,4 +72,19 @@ Foam::scalar Foam::det(SquareMatrix<Type>& matrix)
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
template<class AnyType>
void Foam::SquareMatrix<Type>::operator=(const Identity<AnyType>)
{
Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
for (label i=0; i < this->n(); ++i)
{
this->operator()(i, i) = pTraits<Type>::one;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -52,7 +52,7 @@ template<class Type> class RectangularMatrix;
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
Class SquareMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
@ -65,42 +65,57 @@ public:
// Constructors
//- Null constructor.
//- Construct null
inline SquareMatrix();
//- Construct given number of rows/columns.
inline SquareMatrix(const label n);
//- Construct for given size (rows == cols)
inline explicit SquareMatrix(const label n);
//- Construct from a block of another matrix
template<class MatrixType>
inline SquareMatrix(const ConstMatrixBlock<MatrixType>&);
//- Construct from a block of another matrix
template<class MatrixType>
inline SquareMatrix(const MatrixBlock<MatrixType>&);
//- Construct given number of rows/columns
//- Construct for given size (rows == cols)
//- 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
inline SquareMatrix(const label m, const label n, const zero);
//- Construct given number of rows/columns
//- initializing to the identity matrix
inline SquareMatrix(const label n, const Identity<Type>);
//- Construct with given number of rows and rows
//- Construct for given size (rows == cols)
//- initializing all elements to the given value
inline SquareMatrix(const label n, const Type& val);
//- Construct for given size (rows == cols)
//- initializing to the identity matrix
template<class AnyType>
inline SquareMatrix(const label n, const Identity<AnyType>);
//- Construct given number of rows/columns (checked to be equal)
// For constructor consistency in rectangular matrices
inline explicit SquareMatrix(const labelPair& dims);
//- Construct given number of rows/columns (checked to be equal)
//- initializing all elements to zero
// For constructor consistency with rectangular matrices
inline SquareMatrix(const labelPair& dims, const zero);
//- Construct given number of rows/columns (checked to be equal)
//- initializing all elements to the given value
// For constructor consistency with rectangular matrices
inline SquareMatrix(const labelPair& dims, const Type& val);
//- Construct given number of rows/columns (checked to be equal)
//- initializing all elements to zero
inline SquareMatrix(const label m, const label n, const zero);
//- Construct from sub-matrix block
template<class MatrixType>
inline SquareMatrix(const ConstMatrixBlock<MatrixType>& mat);
//- Construct from sub-matrix block
template<class MatrixType>
inline SquareMatrix(const MatrixBlock<MatrixType>& mat);
//- Construct as copy of a RectangularMatrix
//- which is checked to be square
inline explicit SquareMatrix(const RectangularMatrix<Type>&);
inline explicit SquareMatrix(const RectangularMatrix<Type>& mat);
//- Construct from Istream.
inline SquareMatrix(Istream&);
//- Construct from Istream
inline explicit SquareMatrix(Istream& is);
//- Clone
inline autoPtr<SquareMatrix<Type>> clone() const;
@ -120,16 +135,17 @@ public:
inline void shallowResize(const label m);
// Member operators
// Member Operators
//- Assignment of all elements to zero
void operator=(const zero);
//- Assign all elements to zero
inline void operator=(const zero);
//- Assignment of all elements to the given value
void operator=(const Type& val);
//- Assign all elements to value
inline void operator=(const Type& val);
//- Clear and assign diagonal to identity
void operator=(const Identity<Type>);
//- Set to identity matrix
template<class AnyType>
void operator=(const Identity<AnyType>);
};

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

@ -25,6 +25,15 @@ License
\*---------------------------------------------------------------------------*/
#define CHECK_MATRIX_IS_SQUARE(a, b) \
if (a != b) \
{ \
FatalErrorInFunction \
<< "Attempt to create a non-square matrix (" \
<< a << ", " << b << ')' << nl << abort(FatalError); \
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
@ -41,28 +50,6 @@ inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
{}
template<class Type>
template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const ConstMatrixBlock<MatrixType>& block
)
:
Matrix<SquareMatrix<Type>, Type>(block)
{}
template<class Type>
template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const MatrixBlock<MatrixType>& block
)
:
Matrix<SquareMatrix<Type>, Type>(block)
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
@ -74,26 +61,6 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
{}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label m,
const label n,
const zero
)
:
Matrix<SquareMatrix<Type>, Type>(m, n, Zero)
{
if (m != n)
{
FatalErrorInFunction
<< "Attempt to construct a square matrix "
<< m << " x " << n << nl
<< abort(FatalError);
}
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
@ -106,21 +73,100 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
template<class Type>
template<class AnyType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label n,
const Identity<Type>
const Identity<AnyType>
)
:
Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
{
for (label i=0; i<n; i++)
for (label i=0; i < n; ++i)
{
this->operator()(i, i) = Type(I);
this->operator()(i, i) = pTraits<Type>::one;
}
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const labelPair& dims
)
:
Matrix<SquareMatrix<Type>, Type>(dims)
{
CHECK_MATRIX_IS_SQUARE(dims.first(), dims.second());
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const labelPair& dims,
const zero
)
:
Matrix<SquareMatrix<Type>, Type>(dims, Zero)
{
CHECK_MATRIX_IS_SQUARE(dims.first(), dims.second());
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const labelPair& dims,
const Type& val
)
:
Matrix<SquareMatrix<Type>, Type>(dims, val)
{
CHECK_MATRIX_IS_SQUARE(dims.first(), dims.second());
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label m,
const label n,
const zero
)
:
Matrix<SquareMatrix<Type>, Type>(m, n, Zero)
{
CHECK_MATRIX_IS_SQUARE(m, n);
}
template<class Type>
template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const ConstMatrixBlock<MatrixType>& mat
)
:
Matrix<SquareMatrix<Type>, Type>(mat)
{
// Check is square?
}
template<class Type>
template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const MatrixBlock<MatrixType>& mat
)
:
Matrix<SquareMatrix<Type>, Type>(mat)
{
// Check is square?
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
@ -129,13 +175,7 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
:
Matrix<SquareMatrix<Type>, Type>(mat)
{
if (mat.m() != mat.n())
{
FatalErrorInFunction
<< "Attempt to construct a square matrix from a rectangular matrix "
<< mat.m() << " x " << mat.n() << nl
<< abort(FatalError);
}
CHECK_MATRIX_IS_SQUARE(mat.m(), mat.n());
}
@ -143,7 +183,9 @@ template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
:
Matrix<SquareMatrix<Type>, Type>(is)
{}
{
CHECK_MATRIX_IS_SQUARE(this->m(), this->n());
}
template<class Type>
@ -180,31 +222,19 @@ inline void Foam::SquareMatrix<Type>::shallowResize(const label m)
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
void Foam::SquareMatrix<Type>::operator=(const zero)
inline void Foam::SquareMatrix<Type>::operator=(const zero)
{
Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
}
template<class Type>
void Foam::SquareMatrix<Type>::operator=(const Type& val)
inline 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);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -237,4 +267,9 @@ inline Foam::SquareMatrix<Type> symmOuter
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#undef CHECK_MATRIX_IS_SQUARE
// ************************************************************************* //

17
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.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
@ -110,4 +110,19 @@ Type Foam::det(const SymmetricSquareMatrix<Type>& matrix)
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
template<class AnyType>
void Foam::SymmetricSquareMatrix<Type>::operator=(const Identity<AnyType>)
{
Matrix<SymmetricSquareMatrix<Type>, Type>::operator=(Zero);
for (label i=0; i < this->n(); ++i)
{
this->operator()(i, i) = pTraits<Type>::one;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -40,7 +40,6 @@ SourceFiles
#define SymmetricSquareMatrix_H
#include "SquareMatrix.H"
#include "Identity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -61,27 +60,43 @@ public:
// Constructors
//- Null constructor.
//- Construct null
inline SymmetricSquareMatrix();
//- Construct given number of rows/columns.
inline SymmetricSquareMatrix(const label n);
//- Construct for given size (rows == cols)
inline explicit SymmetricSquareMatrix(const label n);
//- Construct given number of rows/columns, initializing to zero
//- Construct for given size (rows == cols)
//- initializing all elements to zero
inline SymmetricSquareMatrix(const label n, const zero);
//- Construct with given number of rows/columns
//- Construct for given size (rows == cols)
//- 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 for given size (rows == cols)
//- setting the diagonal to one and off-diagonals to zero
template<class AnyType>
inline SymmetricSquareMatrix(const label n, const Identity<AnyType>);
//- Construct from Istream.
inline SymmetricSquareMatrix(Istream&);
//- Construct from Istream
inline explicit SymmetricSquareMatrix(Istream& is);
//- Clone
inline autoPtr<SymmetricSquareMatrix<Type>> clone() const;
// Member Operators
//- Assign all elements to zero
inline void operator=(const zero);
//- Assign all elements to value
inline void operator=(const Type& val);
//- Set to identity matrix
template<class AnyType>
void operator=(const Identity<AnyType>);
};

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

@ -25,6 +25,15 @@ License
\*---------------------------------------------------------------------------*/
#define CHECK_MATRIX_IS_SQUARE(a, b) \
if (a != b) \
{ \
FatalErrorInFunction \
<< "Attempt to create a non-square matrix (" \
<< a << ", " << b << ')' << nl << abort(FatalError); \
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
@ -64,15 +73,16 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
template<class Type>
template<class AnyType>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label n,
const Identity<Type>
const Identity<AnyType>
)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
{
for (label i=0; i<n; i++)
for (label i=0; i < n; ++i)
{
this->operator()(i, i) = pTraits<Type>::one;
}
@ -83,7 +93,9 @@ template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(Istream& is)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(is)
{}
{
CHECK_MATRIX_IS_SQUARE(this->m(), this->n());
}
template<class Type>
@ -94,4 +106,24 @@ Foam::SymmetricSquareMatrix<Type>::clone() const
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
inline void Foam::SymmetricSquareMatrix<Type>::operator=(const zero)
{
Matrix<SymmetricSquareMatrix<Type>, Type>::operator=(Zero);
}
template<class Type>
inline void Foam::SymmetricSquareMatrix<Type>::operator=(const Type& val)
{
Matrix<SymmetricSquareMatrix<Type>, Type>::operator=(val);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#undef CHECK_MATRIX_IS_SQUARE
// ************************************************************************* //

1
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
View File

@ -323,6 +323,7 @@ void Foam::multiply
}
//- Pseudo-inverse algorithm for scalar matrices, using Moore-Penrose inverse
Foam::scalarRectangularMatrix Foam::SVDinv
(
const scalarRectangularMatrix& A,

16
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
View File

@ -60,12 +60,12 @@ typedef SymmetricSquareMatrix<scalar> scalarSymmetricSquareMatrix;
typedef DiagonalMatrix<scalar> scalarDiagonalMatrix;
//- Solve the matrix using Gaussian elimination with pivoting,
// returning the solution in the source
//- returning the solution in the source
template<class Type>
void solve(scalarSquareMatrix& matrix, List<Type>& source);
//- Solve the matrix using Gaussian elimination with pivoting
// and return the solution
//- and return the solution
template<class Type>
void solve
(
@ -82,7 +82,7 @@ void LUDecompose
);
//- LU decompose the matrix with pivoting.
// sign is -1 for odd number of row interchanges and 1 for even number.
//- sign is -1 for odd number of row interchanges and 1 for even number.
void LUDecompose
(
scalarSquareMatrix& matrix,
@ -94,7 +94,7 @@ void LUDecompose
void LUDecompose(scalarSymmetricSquareMatrix& matrix);
//- LU back-substitution with given source, returning the solution
// in the source
//- in the source
template<class Type>
void LUBacksubstitute
(
@ -104,8 +104,8 @@ void LUBacksubstitute
);
//- LU back-substitution with given source, returning the solution
// in the source. Specialised for symmetric square matrices that have been
// decomposed into LU where U = L.T() as pivoting is not required
//- in the source. Specialised for symmetric square matrices that have been
//- decomposed into LU where U = L.T() as pivoting is not required
template<class Type>
void LUBacksubstitute
(
@ -114,12 +114,12 @@ void LUBacksubstitute
);
//- Solve the matrix using LU decomposition with pivoting
// returning the LU form of the matrix and the solution in the source
//- returning the LU form of the matrix and the solution in the source
template<class Type>
void LUsolve(scalarSquareMatrix& matrix, List<Type>& source);
//- Solve the matrix using LU decomposition returning the LU form of the matrix
// and the solution in the source, where U = L.T()
//- and the solution in the source, where U = L.T()
template<class Type>
void LUsolve(scalarSymmetricSquareMatrix& matrix, List<Type>& source);