diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index ad48481fa1..95cfa4c1af 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -25,6 +25,7 @@ License
 
 #include "scalarMatrices.H"
 #include "vector.H"
+#include "tensor.H"
 #include "IFstream.H"
 
 using namespace Foam;
@@ -50,7 +51,7 @@ int main(int argc, char *argv[])
     Info<< max(hmm) << endl;
     Info<< min(hmm) << endl;
 
-    SquareMatrix<scalar> hmm2(3, 3, 1.0);
+    SquareMatrix<scalar> hmm2(3, I);
 
     hmm = hmm2;
 
@@ -72,7 +73,14 @@ int main(int argc, char *argv[])
     Info<< hmm5 << endl;
 
     {
-        scalarSymmetricSquareMatrix symmMatrix(3, 3, 0);
+        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;
+    }
+
+    {
+        scalarSymmetricSquareMatrix symmMatrix(3, Zero);
 
         symmMatrix(0, 0) = 4;
         symmMatrix(1, 0) = 12;
@@ -104,7 +112,7 @@ int main(int argc, char *argv[])
     }
 
     {
-        scalarSquareMatrix squareMatrix(3, 3, 0);
+        scalarSquareMatrix squareMatrix(3, Zero);
 
         squareMatrix(0, 0) = 4;
         squareMatrix(0, 1) = 12;
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index ed00d1a6cf..d5f295a62b 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -675,7 +675,6 @@ int main(int argc, char *argv[])
     // Matrix sum in j(Fij) for each i (if enclosure sum = 1)
     scalarSquareMatrix sumViewFactorPatch
     (
-        totalPatches,
         totalPatches,
         0.0
     );
@@ -888,7 +887,7 @@ int main(int argc, char *argv[])
 
     if (Pstream::master())
     {
-        scalarSquareMatrix Fmatrix(totalNCoarseFaces, totalNCoarseFaces, 0.0);
+        scalarSquareMatrix Fmatrix(totalNCoarseFaces, 0.0);
 
         labelListList globalFaceFaces(visibleFaceFaces.size());
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 510305f28c..2882e6af2a 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -303,7 +303,7 @@ triSurfacePointScalarField calcCurvature
         faceCoordSys.normalize();
 
         // Construct the matrix to solve
-        scalarSymmetricSquareMatrix T(3, 3, 0);
+        scalarSymmetricSquareMatrix T(3, 0);
         scalarDiagonalMatrix Z(3, 0);
 
         // Least Squares
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index 36acb8a6b1..a095a7c507 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,7 +50,7 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     a_(iMaxX_, 0.0),
-    alpha_(kMaxX_, kMaxX_, 0.0),
+    alpha_(kMaxX_, 0.0),
     d_p_(n_, kMaxX_, 0.0),
     x_p_(kMaxX_, 0.0),
     err_(kMaxX_, 0.0),
@@ -60,7 +60,7 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
     yErr_(n_, 0.0),
     dydx0_(n_),
     dfdx_(n_, 0.0),
-    dfdy_(n_, n_, 0.0),
+    dfdy_(n_, 0.0),
     first_(1),
     epsOld_(-1.0)
 {}
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index 74439c15de..cf34d78a21 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
@@ -35,7 +35,7 @@ Foam::scalar Foam::detDecomposed
     const label sign
 )
 {
-    scalar diagProduct = 1.0;
+    Type diagProduct = pTraits<Type>::one;
 
     for (label i = 0; i < matrix.m(); ++i)
     {
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index 270891a955..39403194d5 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -38,6 +38,7 @@ SourceFiles
 #define SquareMatrix_H
 
 #include "Matrix.H"
+#include "Identity.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,14 +65,17 @@ public:
         //- Construct given number of rows/columns.
         inline SquareMatrix(const label n);
 
-        //- Construct given number of rows and columns,
-        //  It checks that m == n.
-        inline SquareMatrix(const label m, const label n);
+        //- Construct given number of rows/columns
+        //  initializing all elements to zero
+        inline SquareMatrix(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
-        //  and value for all elements.
-        //  It checks that m == n.
-        inline SquareMatrix(const label m, const label n, const Type&);
+        //  initializing all elements to the given value
+        inline SquareMatrix(const label n, const Type&);
 
         //- Construct from Istream.
         inline SquareMatrix(Istream&);
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 74d7bdb95c..ffa595d51d 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -31,40 +31,51 @@ inline Foam::SquareMatrix<Type>::SquareMatrix()
     Matrix<SquareMatrix<Type>, Type>()
 {}
 
+
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
 :
     Matrix<SquareMatrix<Type>, Type>(n, n)
 {}
 
-template<class Type>
-inline Foam::SquareMatrix<Type>::SquareMatrix(const label m, const label n)
-:
-    Matrix<SquareMatrix<Type>, Type>(m, n)
-{
-    if (m != n)
-    {
-        FatalErrorInFunction
-          << "m != n for constructing a square matrix" << exit(FatalError);
-    }
-}
 
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix
 (
-    const label m,
+    const label n,
+    const zero
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
+{}
+
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+    const label n,
+    const Identity<Type>
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
+{
+    for (label i = 0; i < n; ++i)
+    {
+        this->operator()(i, i) = I;
+    }
+}
+
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
     const label n,
     const Type& t
 )
 :
-    Matrix<SquareMatrix<Type>, Type>(m, n, t)
-{
-    if (m != n)
-    {
-        FatalErrorInFunction
-          << "m != n for constructing a square matrix" << exit(FatalError);
-    }
-}
+    Matrix<SquareMatrix<Type>, Type>(n, n, t)
+{}
+
 
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
@@ -72,6 +83,7 @@ inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
     Matrix<SquareMatrix<Type>, Type>(is)
 {}
 
+
 template<class Type>
 inline Foam::autoPtr<Foam::SquareMatrix<Type>>
 Foam::SquareMatrix<Type>::clone() const
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
index 28ad1d2745..dcb31ef241 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
@@ -33,15 +33,17 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
     const SymmetricSquareMatrix<Type>& matrix
 )
 {
-    SymmetricSquareMatrix<Type> inv(matrix.m(), matrix.m(), 0.0);
+    const label n = matrix.n();
 
-    for (label i = 0; i < matrix.m(); ++i)
+    SymmetricSquareMatrix<Type> inv(n, Zero);
+
+    for (label i = 0; i < n; ++i)
     {
         inv(i, i) = 1.0/matrix(i, i);
 
         for (label j = 0; j < i; ++j)
         {
-            scalar sum = 0.0;
+            Type sum = Zero;
 
             for (label k = j; k < i; k++)
             {
@@ -52,7 +54,20 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
         }
     }
 
-    return inv.T()*inv;
+    SymmetricSquareMatrix<Type> result(n, Zero);
+
+    for (label k = 0; k < n; k++)
+    {
+        for (label i = 0; i <= k; i++)
+        {
+            for (label j = 0; j <= k; j++)
+            {
+                result(i, j) += inv(k, i)*inv(k, j);
+            }
+        }
+    }
+
+    return result;
 }
 
 
@@ -63,7 +78,6 @@ Foam::SymmetricSquareMatrix<Type> Foam::inv
 )
 {
     SymmetricSquareMatrix<Type> matrixTmp(matrix);
-
     LUDecompose(matrixTmp);
 
     return invDecomposed(matrixTmp);
@@ -71,9 +85,9 @@ Foam::SymmetricSquareMatrix<Type> Foam::inv
 
 
 template<class Type>
-Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
+Type Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
 {
-    scalar diagProduct = 1.0;
+    Type diagProduct = pTraits<Type>::one;
 
     for (label i = 0; i < matrix.m(); ++i)
     {
@@ -85,10 +99,9 @@ Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
 
 
 template<class Type>
-Foam::scalar Foam::det(const SymmetricSquareMatrix<Type>& matrix)
+Type Foam::det(const SymmetricSquareMatrix<Type>& matrix)
 {
-    SymmetricSquareMatrix<Type> matrixTmp = matrix;
-
+    SymmetricSquareMatrix<Type> matrixTmp(matrix);
     LUDecompose(matrixTmp);
 
     return detDecomposed(matrixTmp);
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
index 193ee1602f..543d16b848 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
@@ -38,6 +38,7 @@ SourceFiles
 #define SymmetricSquareMatrix_H
 
 #include "SquareMatrix.H"
+#include "Identity.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,12 +65,15 @@ public:
         //- Construct given number of rows/columns.
         inline SymmetricSquareMatrix(const label n);
 
-        //- Construct with given number of rows/columns
-        inline SymmetricSquareMatrix(const label m, const label n);
+        //- Construct given number of rows/columns, initializing to zero
+        inline SymmetricSquareMatrix(const label n, const zero);
+
+        //- Construct given number of rows/columns,
+        inline SymmetricSquareMatrix(const label n, const Identity<Type>);
 
         //- Construct with given number of rows/columns
-        //  and value for all elements.
-        inline SymmetricSquareMatrix(const label m, const label n, const Type&);
+        //  initializing all elements to the given value
+        inline SymmetricSquareMatrix(const label n, const Type&);
 
         //- Construct from Istream.
         inline SymmetricSquareMatrix(Istream&);
@@ -91,11 +95,11 @@ SymmetricSquareMatrix<Type> inv(const SymmetricSquareMatrix<Type>&);
 
 //- Return the LU decomposed SymmetricSquareMatrix det
 template<class Type>
-scalar detDecomposed(const SymmetricSquareMatrix<Type>&);
+Type detDecomposed(const SymmetricSquareMatrix<Type>&);
 
 //- Return the SymmetricSquareMatrix det
 template<class Type>
-scalar det(const SymmetricSquareMatrix<Type>&);
+Type det(const SymmetricSquareMatrix<Type>&);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
index a363026aac..821daf199b 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
@@ -42,17 +42,26 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(const label n)
 template<class Type>
 inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
 (
-    const label m,
-    const label n
+    const label n,
+    const zero
 )
 :
-    Matrix<SymmetricSquareMatrix<Type>, Type>(m, n)
+    Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
+{}
+
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
+(
+    const label n,
+    const Identity<Type>
+)
+:
+    Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
 {
-    if (m != n)
+    for (label i = 0; i < n; i++)
     {
-        FatalErrorInFunction
-            << "m != n for constructing a symmetric square matrix"
-            << exit(FatalError);
+        this->operator()(i, i) = pTraits<Type>::one;
     }
 }
 
@@ -60,20 +69,12 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
 template<class Type>
 inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
 (
-    const label m,
     const label n,
     const Type& t
 )
 :
-    Matrix<SymmetricSquareMatrix<Type>, Type>(m, n, t)
-{
-    if (m != n)
-    {
-        FatalErrorInFunction
-            << "m != n for constructing a symmetric square matrix"
-            << exit(FatalError);
-    }
-}
+    Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, t)
+{}
 
 
 template<class Type>
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
index 739251c358..7ac02d189f 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
@@ -35,7 +35,7 @@ Foam::tmp<Foam::Field<Type>> Foam::lduMatrix::H(const Field<Type>& psi) const
 {
     tmp<Field<Type>> tHpsi
     (
-        new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
+        new Field<Type>(lduAddr().size(), Zero)
     );
 
     if (lowerPtr_ || upperPtr_)
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
index c6acc739e2..52413a38dd 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -87,7 +87,7 @@ void Foam::solve
     // Back-substitution
     for (label j=m-1; j>=0; j--)
     {
-        Type ntempvec = pTraits<Type>::zero;
+        Type ntempvec = Zero;
 
         for (label k=j+1; k<m; k++)
         {
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
index 960ba0d563..0f1c6995b9 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
@@ -43,7 +43,7 @@ Foam::simpleMatrix<Type>::simpleMatrix
     const Type& sourceVal
 )
 :
-    scalarSquareMatrix(mSize, mSize, coeffVal),
+    scalarSquareMatrix(mSize, coeffVal),
     source_(mSize, sourceVal)
 {}
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 2c8a7e2ba1..d70eccfa9c 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -173,7 +173,7 @@ void Foam::radiation::viewFactor::initialise()
     {
         Fmatrix_.reset
         (
-            new scalarSquareMatrix(totalNCoarseFaces_, totalNCoarseFaces_, 0.0)
+            new scalarSquareMatrix(totalNCoarseFaces_, 0.0)
         );
 
         if (debug)
@@ -224,12 +224,7 @@ void Foam::radiation::viewFactor::initialise()
         {
             CLU_.reset
             (
-                new scalarSquareMatrix
-                (
-                    totalNCoarseFaces_,
-                    totalNCoarseFaces_,
-                    0.0
-                )
+                new scalarSquareMatrix(totalNCoarseFaces_, 0.0)
             );
 
             pivotIndices_.setSize(CLU_().m());
@@ -529,7 +524,7 @@ void Foam::radiation::viewFactor::calculate()
         // Variable emissivity
         if (!constEmissivity_)
         {
-            scalarSquareMatrix C(totalNCoarseFaces_, totalNCoarseFaces_, 0.0);
+            scalarSquareMatrix C(totalNCoarseFaces_, 0.0);
 
             for (label i=0; i<totalNCoarseFaces_; i++)
             {