diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index a402c73a30..4e001d7c6a 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 0c0c6b1f1b..8c39029e4e 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 f03b4025f2..10b2e41b68 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 ed6fd895b6..1b0908dc3c 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 9de6f4f3f0..6024d21c18 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 e1037eec5d..d6e3d836c6 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 d2b3d119ac..27e90c7a9f 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 9ac2f2977e..408f1393a1 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 5b5e68e8b6..68ecd56a59 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           |
+    \\  /    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);
 };
 
 
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
index 26c28ccfac..da3e03b53f 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           |
+    \\  /    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
 
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index 5429f5c78e..b48d01545c 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 44893b92a6..dec80275bd 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 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>);
 };
 
 
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 0d5a3a8329..4e1b3ca483 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>
@@ -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
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
index 82ef280784..561cd7bafc 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 f4bf7c9c41..7390f401ab 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 a2f8dc0109..c707347a06 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 e6e1ca3e85..b8551c4ee2 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 19f0705ed8..48990578c9 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);