Commit a6e11607 authored by Mark Olesen's avatar Mark Olesen

ENH: harmonize matrix constructors (#1220)

- 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
parent 773d72e0
......@@ -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;
......
......@@ -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);
}
......
......@@ -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
);
......
......@@ -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: ("
......
......@@ -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;
......
......@@ -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
{
......
......@@ -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()
(
......
......@@ -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>
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
......@@ -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 given number of rows/columns
//- initializing all elements to zero
inline RectangularMatrix(const label m, const label n, const zero);
//- 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>&);
inline RectangularMatrix(const ConstMatrixBlock<MatrixType>& mat);
//- Construct from a block of another matrix
template<class MatrixType>
inline RectangularMatrix(const MatrixBlock<MatrixType>&);
//- Construct with given number of rows and 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.
inline RectangularMatrix(const label m, const label n, const Type& val);
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);
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
......@@ -46,24 +46,24 @@ 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)
Matrix<RectangularMatrix<Type>, Type>(n, n)
{}
template<class Type>
template<class MatrixType>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const MatrixBlock<MatrixType>& block
const label m,
const label n,
const zero
)
:
Matrix<RectangularMatrix<Type>, Type>(block)
Matrix<RectangularMatrix<Type>, Type>(m, n, Zero)
{}
......@@ -72,22 +72,64 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const label m,
const label n,
const Type& val
)
:
Matrix<RectangularMatrix<Type>, Type>(m, n, val)
{}
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
)
:
Matrix<RectangularMatrix<Type>, Type>(m, n, Zero)
RectangularMatrix<Type>(dims.first(), dims.second(), Zero)
{}
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
const label m,
const label n,
const labelPair& dims,
const Type& val
)
:
Matrix<RectangularMatrix<Type>, Type>(m, n, 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)
{}
......@@ -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
// ************************************************************************* //
......@@ -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;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -52,7 +52,7 @@ template<class Type> class RectangularMatrix;
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
Class SquareMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
......@@ -65,42 +65,57 @@ public:
// Constructors