diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index 61aaf2d8a35ae0c89c3f9f38a8e7004cf4b183a7..d53bfd2e75c28ace421447e8598d8a1584bd562a 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -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;
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
index 84a6c32efccc258250a93ed4562c9c5016fa591e..f82f2c137f6c770f8e162c77b9d4edf0f2f9d2ab 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
@@ -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);
}
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
index c5fb4c48e46d265cd59481e10ad129dec70ab509..0c65a4cf9f8fc4c80dd679cabcf98055b7612747 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
@@ -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
);
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index ba1393b8807d9be7345f5acd42c27e276850a2c2..204899f33ec777300115c6d6e80582fad924f99c 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -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: ("
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 07dcdeca022bb8ef290d2eefcee430d453c4e14a..a9fa0c3fe9b819dc82a177da8b625e04f772300e 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -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;
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 315ac7770dd37c5c0951548e77af77676b4ee547..40cb96e99969c0a645aba800bceca88864ba4bdf 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -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
{
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
index d2b3d119aca855c215c2de33b72c81aba6d5d5c4..27e90c7a9f4d77811d82e8fcd720a17201734e86 100644
--- a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd |
+ \\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016 OpenFOAM Foundation
@@ -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()
(
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
index 9ac2f2977e8da8d9f1dddb97d8f74a12eede9296..408f1393a10213472dbbc78186f4e129f6b59829 100644
--- a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd |
+ \\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016 OpenFOAM Foundation
@@ -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>
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
index ae49b677e51d57c4eec8edb09493b923e183a8a2..0d1474034c9313d3cfdcdb019968db54491ff0d8 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -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);
};
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
index f45b3b797803ee4625265b241b29c0dcd74dd105..b50f37ef202df443070127383e2290c1e85f22c4 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -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
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index 5429f5c78e85d18b71907ff30299533797ff3e2d..b48d01545c376a421fde9fb38c1d4d9f4afa925f 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
@@ -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;
+ }
+}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index f1819a9c9cb46385ddacb248f52b64ef4589a733..8ff429ed3e1e2186286e123e0586dc24f477b077 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -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 for given size (rows == cols)
+ //- initializing all elements to zero
+ inline SquareMatrix(const label n, const zero);
- //- Construct from a block of another matrix
- template<class MatrixType>
- inline SquareMatrix(const MatrixBlock<MatrixType>&);
+ //- 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
+ //- Construct given number of rows/columns (checked to be equal)
//- initializing all elements to zero
- inline SquareMatrix(const label n, const 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 and columns (checked to be equal)
+ //- 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 given number of rows/columns
- //- initializing to the identity matrix
- inline SquareMatrix(const label n, const Identity<Type>);
+ //- Construct from sub-matrix block
+ template<class MatrixType>
+ inline SquareMatrix(const ConstMatrixBlock<MatrixType>& mat);
- //- Construct with given number of rows and rows
- //- initializing all elements to the given value
- inline SquareMatrix(const label n, const Type& val);
+ //- 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>);
};
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 97c61c892372a0da020ec6d1010dabb497572ec7..9cba19c434ced0f93055ea0dbcff776ced27a5c2 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -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>
@@ -42,36 +51,80 @@ inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
template<class Type>
-template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
- const ConstMatrixBlock<MatrixType>& block
+ const label n,
+ const zero
)
:
- Matrix<SquareMatrix<Type>, Type>(block)
+ Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
{}
template<class Type>
-template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
- const MatrixBlock<MatrixType>& block
+ const label n,
+ const Type& val
)
:
- Matrix<SquareMatrix<Type>, Type>(block)
+ Matrix<SquareMatrix<Type>, Type>(n, n, val)
{}
template<class Type>
+template<class AnyType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const label n,
- const zero
+ const Identity<AnyType>
)
:
Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
-{}
+{
+ for (label i=0; i < n; ++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>
@@ -84,40 +137,33 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
:
Matrix<SquareMatrix<Type>, Type>(m, n, Zero)
{
- if (m != n)
- {
- FatalErrorInFunction
- << "Attempt to construct a square matrix "
- << m << " x " << n << nl
- << abort(FatalError);
- }
+ CHECK_MATRIX_IS_SQUARE(m, n);
}
template<class Type>
+template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
- const label n,
- const Type& val
+ const ConstMatrixBlock<MatrixType>& mat
)
:
- Matrix<SquareMatrix<Type>, Type>(n, n, val)
-{}
+ Matrix<SquareMatrix<Type>, Type>(mat)
+{
+ // Check is square?
+}
template<class Type>
+template<class MatrixType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
- const label n,
- const Identity<Type>
+ const MatrixBlock<MatrixType>& mat
)
:
- Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
+ Matrix<SquareMatrix<Type>, Type>(mat)
{
- for (label i=0; i<n; i++)
- {
- this->operator()(i, i) = Type(I);
- }
+ // Check is square?
}
@@ -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
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
index 82ef2807843b6c3f6e0043c1b1294eb686b12db4..561cd7bafc0bce3c4bea49da4c51c007b2fc45f0 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd |
+ \\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 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;
+ }
+}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
index f4bf7c9c41aaf24574b79b8528877e8f464d9bd7..7390f401ab54f9f5cd2810b7a2c8cdeedbefda24 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
@@ -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>);
};
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
index a2f8dc01092d2bc367084a4f1f90c581350faedf..c707347a06c11b42af217884272429c6235b3ced 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
@@ -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
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index 9b7533e6183fb04bb8a2ca2d4a3fbe89e2528b25..c473c9f490fbb263cedbe0df09bcc2e2c6a70b10 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
@@ -323,6 +323,7 @@ void Foam::multiply
}
+//- Pseudo-inverse algorithm for scalar matrices, using Moore-Penrose inverse
Foam::scalarRectangularMatrix Foam::SVDinv
(
const scalarRectangularMatrix& A,
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
index a8b026a50519fd6f5cdfc9f2d48650758cc2dd75..8c01c02a07a0d26105e997cefbe8a2b83fe0e332 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
@@ -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);