ENH: partial overhaul of Matrix type (#1220)

- additional operators:
  + compound assignment
  + inner product: operator&
  + outer product: operator^

- additional functions:
   - MatrixBlock methods: subColumn, subRow, subMatrix
   - L2 norms for matrix or column
   - trace, diag, round, transpose

- MatrixBlock methods: col(), block() are deprecated since their
  access patterns with (size, offset) are unnatural/unwieldy.

- verifications by test/Matrix/Test-Matrix

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -334,16 +334,29 @@ void Foam::Matrix<Form, Type>::resize(const label m, const label n)
 }
 
 
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::round(const scalar tol)
+{
+    for (Type& val : *this)
+    {
+        if (mag(val) < tol)
+        {
+            val = Zero;
+        }
+    }
+}
+
+
 template<class Form, class Type>
 Form Foam::Matrix<Form, Type>::T() const
 {
-    Form At(n(), m());
+    Form At(labelPair{n(), m()});
 
     for (label i = 0; i < m(); ++i)
     {
         for (label j = 0; j < n(); ++j)
         {
-            At(j, i) = (*this)(i, j);
+            At(j, i) = Detail::conj((*this)(i, j));
         }
     }
 
@@ -351,6 +364,91 @@ Form Foam::Matrix<Form, Type>::T() const
 }
 
 
+template<class Form, class Type>
+Foam::List<Type> Foam::Matrix<Form, Type>::diag() const
+{
+    const label len = Foam::min(mRows_, nCols_);
+
+    List<Type> result(len);
+
+    for (label i=0; i < len; ++i)
+    {
+        result[i] = (*this)(i, i);
+    }
+
+    return result;
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::diag(const UList<Type>& list)
+{
+    const label len = Foam::min(mRows_, nCols_);
+
+    #ifdef FULLDEBUG
+    if (list.size() != len)
+    {
+        FatalErrorInFunction
+            << "List size (" << list.size()
+            << ") incompatible with Matrix diagonal" << abort(FatalError);
+    }
+    #endif
+
+    for (label i=0; i < len; ++i)
+    {
+        (*this)(i, i) = list[i];
+    }
+}
+
+
+template<class Form, class Type>
+Type Foam::Matrix<Form, Type>::trace() const
+{
+    const label len = Foam::min(mRows_, nCols_);
+
+    Type val = Zero;
+
+    for (label i=0; i < len; ++i)
+    {
+        val += (*this)(i, i);
+    }
+
+    return val;
+}
+
+
+template<class Form, class Type>
+Foam::scalar Foam::Matrix<Form, Type>::columnNorm
+(
+    const label colIndex,
+    const bool noSqrt
+) const
+{
+    scalar result = Zero;
+
+    for (label i=0; i < mRows_; ++i)
+    {
+        result += magSqr((*this)(i, colIndex));
+    }
+
+    return noSqrt ? result : Foam::sqrt(result);
+}
+
+
+template<class Form, class Type>
+Foam::scalar Foam::Matrix<Form, Type>::norm(const bool noSqrt) const
+{
+    scalar result = Zero;
+
+    for (const Type& val : *this)
+    {
+        result += magSqr(val);
+    }
+
+    return noSqrt ? result : Foam::sqrt(result);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Form, class Type>
@@ -416,9 +514,9 @@ void Foam::Matrix<Form, Type>::operator=
     const MatrixBlock<MatrixType>& Mb
 )
 {
-    for (label i=0; i < mRows_; ++i)
+    for (label i = 0; i < mRows_; ++i)
     {
-        for (label j=0; j < nCols_; ++j)
+        for (label j = 0; j < nCols_; ++j)
         {
             (*this)(i, j) = Mb(i, j);
         }
@@ -443,6 +541,7 @@ void Foam::Matrix<Form, Type>::operator=(const zero)
 template<class Form, class Type>
 void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& other)
 {
+    #ifdef FULLDEBUG
     if (this == &other)
     {
         FatalErrorInFunction
@@ -458,15 +557,13 @@ void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& other)
             << other.m() << ", " << other.n() << ')' << nl
             << abort(FatalError);
     }
+    #endif
 
-    Type* out = this->data();
+    auto iter2 = other.cbegin();
-    const Type* in = other.cdata();
+    for (Type& val : *this)
-
-    const label len = this->size();
-
-    for (label idx = 0; idx < len; ++idx)
     {
-        out[idx] += in[idx];
+        val += *iter2;
+        ++iter2;
     }
 }
 
@@ -474,6 +571,7 @@ void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& other)
 template<class Form, class Type>
 void Foam::Matrix<Form, Type>::operator-=(const Matrix<Form, Type>& other)
 {
+    #ifdef FULLDEBUG
     if (this == &other)
     {
         FatalErrorInFunction
@@ -489,21 +587,39 @@ void Foam::Matrix<Form, Type>::operator-=(const Matrix<Form, Type>& other)
             << other.m() << ", " << other.n() << ')' << nl
             << abort(FatalError);
     }
+    #endif
 
-    Type* out = this->data();
+    auto iter2 = other.cbegin();
-    const Type* in = other.cdata();
+    for (Type& val : *this)
-
-    const label len = this->size();
-
-    for (label idx=0; idx < len; ++idx)
     {
-        out[idx] -= in[idx];
+        val -= *iter2;
+        ++iter2;
     }
 }
 
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator*=(const scalar s)
+void Foam::Matrix<Form, Type>::operator+=(const Type& s)
+{
+    for (Type& val : *this)
+    {
+        val += s;
+    }
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator-=(const Type& s)
+{
+    for (Type& val : *this)
+    {
+        val -= s;
+    }
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator*=(const Type& s)
 {
     for (Type& val : *this)
     {
@@ -513,7 +629,7 @@ void Foam::Matrix<Form, Type>::operator*=(const scalar s)
 
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator/=(const scalar s)
+void Foam::Matrix<Form, Type>::operator/=(const Type& s)
 {
     for (Type& val : *this)
     {
@@ -522,7 +638,7 @@ void Foam::Matrix<Form, Type>::operator/=(const scalar s)
 }
 
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
 
 template<class Form, class Type>
 const Type& Foam::max(const Matrix<Form, Type>& mat)
@@ -569,7 +685,7 @@ Foam::MinMax<Type> Foam::minMax(const Matrix<Form, Type>& mat)
 template<class Form, class Type>
 Form Foam::operator-(const Matrix<Form, Type>& mat)
 {
-    Form result(mat.m(), mat.n());
+    Form result(mat.sizes());
 
     std::transform
     (
@@ -583,9 +699,14 @@ Form Foam::operator-(const Matrix<Form, Type>& mat)
 }
 
 
-template<class Form, class Type>
+template<class Form1, class Form2, class Type>
-Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
+Form1 Foam::operator+
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& B
+)
 {
+    #ifdef FULLDEBUG
     if (A.m() != B.m() || A.n() != B.n())
     {
         FatalErrorInFunction
@@ -594,27 +715,31 @@ Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
             << B.m() << ", " << B.n() << ')' << nl
             << abort(FatalError);
     }
+    #endif
 
-    Form AB(A.m(), A.n());
+    Form1 result(A.sizes());
 
-    Type* ABv = AB.data();
+    std::transform
-    const Type* Av = A.cdata();
+    (
-    const Type* Bv = B.cdata();
+        A.cbegin(),
+        A.cend(),
+        B.cbegin(),
+        result.begin(),
+        std::plus<Type>()
+    );
 
-    const label len = A.size();
+    return result;
-
-    for (label idx = 0; idx < len; ++idx)
-    {
-        ABv[idx] = Av[idx] + Bv[idx];
-    }
-
-    return AB;
 }
 
 
-template<class Form, class Type>
+template<class Form1, class Form2, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
+Form1 Foam::operator-
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& B
+)
 {
+    #ifdef FULLDEBUG
     if (A.m() != B.m() || A.n() != B.n())
     {
         FatalErrorInFunction
@@ -623,85 +748,117 @@ Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
             << B.m() << ", " << B.n() << ')' << nl
             << abort(FatalError);
     }
+    #endif
 
-    Form AB(A.m(), A.n());
+    Form1 result(A.sizes());
 
-    Type* ABv = AB.data();
+    std::transform
-    const Type* Av = A.cdata();
+    (
-    const Type* Bv = B.cdata();
+        A.cbegin(),
-
+        A.cend(),
-    const label len = A.size();
+        B.cbegin(),
-
+        result.begin(),
-    for (label idx=0; idx < len; ++idx)
+        std::minus<Type>()
-    {
+    );
-        ABv[idx] = Av[idx] - Bv[idx];
-    }
-
-    return AB;
-}
-
-
-template<class Form, class Type>
-Form Foam::operator*(const scalar s, const Matrix<Form, Type>& mat)
-{
-    Form result(mat.m(), mat.n());
-
-    const label len = mat.size();
-
-    if (len)
-    {
-        Type* out = result.data();
-        const Type* in = mat.cdata();
-
-        for (label idx = 0; idx < len; ++idx)
-        {
-            out[idx] = s * in[idx];
-        }
-    }
 
     return result;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& mat, const scalar s)
+Form Foam::operator*(const Type& s, const Matrix<Form, Type>& mat)
 {
-    Form result(mat.m(), mat.n());
+    Form result(mat.sizes());
 
-    const label len = mat.size();
+    std::transform
-
+    (
-    if (len)
+        mat.cbegin(),
-    {
+        mat.cend(),
-        Type* out = result.data();
+        result.begin(),
-        const Type* in = mat.cdata();
+        [&](const Type& val) { return s * val; }
-
+    );
-        for (label idx=0; idx < len; ++idx)
-        {
-            out[idx] = in[idx] * s;
-        }
-    }
 
     return result;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator/(const Matrix<Form, Type>& mat, const scalar s)
+Form Foam::operator+(const Type& s, const Matrix<Form, Type>& mat)
 {
-    Form result(mat.m(), mat.n());
+    Form result(mat.sizes());
 
-    const label len = mat.size();
+    std::transform
+    (
+        mat.cbegin(),
+        mat.cend(),
+        result.begin(),
+        [&](const Type& val) { return s + val; }
+    );
 
-    if (len)
+    return result;
-    {
+}
-        Type* out = result.data();
-        const Type* in = mat.cdata();
 
-        for (label idx=0; idx < len; ++idx)
+
-        {
+template<class Form, class Type>
-            out[idx] = in[idx] / s;
+Form Foam::operator-(const Type& s, const Matrix<Form, Type>& mat)
-        }
+{
-    }
+    Form result(mat.sizes());
+
+    std::transform
+    (
+        mat.cbegin(),
+        mat.end(),
+        result.begin(),
+        [&](const Type& val) { return s - val; }
+    );
+
+    return result;
+}
+
+
+template<class Form, class Type>
+Form Foam::operator*(const Matrix<Form, Type>& mat, const Type& s)
+{
+    return s*mat;
+}
+
+
+template<class Form, class Type>
+Form Foam::operator+(const Matrix<Form, Type>& mat, const Type& s)
+{
+    return s + mat;
+}
+
+
+template<class Form, class Type>
+Form Foam::operator-(const Matrix<Form, Type>& mat, const Type& s)
+{
+    Form result(mat.sizes());
+
+    std::transform
+    (
+        mat.cbegin(),
+        mat.end(),
+        result.begin(),
+        [&](const Type& val) { return val - s; }
+    );
+
+    return result;
+}
+
+
+template<class Form, class Type>
+Form Foam::operator/(const Matrix<Form, Type>& mat, const Type& s)
+{
+    Form result(mat.sizes());
+
+    std::transform
+    (
+        mat.cbegin(),
+        mat.end(),
+        result.begin(),
+        [&](const Type& val) { return val / s; }
+    );
 
     return result;
 }
@@ -715,6 +872,7 @@ Foam::operator*
     const Matrix<Form2, Type>& B
 )
 {
+    #ifdef FULLDEBUG
     if (A.n() != B.m())
     {
         FatalErrorInFunction
@@ -724,6 +882,7 @@ Foam::operator*
             << "The columns of A must equal rows of B"
             << abort(FatalError);
     }
+    #endif
 
     typename typeOfInnerProduct<Type, Form1, Form2>::type AB
     (
@@ -732,13 +891,97 @@ Foam::operator*
         Zero
     );
 
-    for (label i=0; i < AB.m(); ++i)
+    for (label i = 0; i < AB.m(); ++i)
     {
-        for (label j=0; j < AB.n(); ++j)
+        for (label k = 0; k < B.m(); ++k)
         {
-            for (label k=0; k < B.m(); ++k)
+            for (label j = 0; j < AB.n(); ++j)
             {
-                AB(i, j) += A(i, k) * B(k, j);
+                AB(i, j) += A(i, k)*B(k, j);
+            }
+        }
+    }
+
+    return AB;
+}
+
+
+template<class Form1, class Form2, class Type>
+typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+Foam::operator&
+(
+    const Matrix<Form1, Type>& AT,
+    const Matrix<Form2, Type>& B
+)
+{
+    #ifdef FULLDEBUG
+    if (AT.m() != B.m())
+    {
+        FatalErrorInFunction
+            << "Attempt to multiply incompatible matrices:" << nl
+            << "Matrix A : (" << AT.m() << ", " << AT.n() << ')' << nl
+            << "Matrix B : (" << B.m() << ", " << B.n() << ')' << nl
+            << "The rows of A must equal rows of B"
+            << abort(FatalError);
+    }
+    #endif
+
+    typename typeOfInnerProduct<Type, Form1, Form2>::type AB
+    (
+        AT.n(),
+        B.n(),
+        Zero
+    );
+
+    for (label i = 0; i < AB.m(); ++i)
+    {
+        for (label k = 0; k < B.m(); ++k)
+        {
+            for (label j = 0; j < AB.n(); ++j)
+            {
+                AB(i, j) += Detail::conj(AT(k, i))*B(k, j);
+            }
+        }
+    }
+
+    return AB;
+}
+
+
+template<class Form1, class Form2, class Type>
+typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+Foam::operator^
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& BT
+)
+{
+    #ifdef FULLDEBUG
+    if (A.n() != BT.n())
+    {
+        FatalErrorInFunction
+            << "Attempt to multiply incompatible matrices:" << nl
+            << "Matrix A : (" << A.m() << ", " << A.n() << ')' << nl
+            << "Matrix B : (" << BT.m() << ", " << BT.n() << ')' << nl
+            << "The columns of A must equal columns of B"
+            << abort(FatalError);
+    }
+    #endif
+
+    typename typeOfInnerProduct<Type, Form1, Form2>::type AB
+    (
+        A.m(),
+        BT.m(),
+        Zero
+    );
+
+    for (label i = 0; i < AB.m(); ++i)
+    {
+        for (label k = 0; k < BT.n(); ++k)
+        {
+            for (label j = 0; j < AB.n(); ++j)
+            {
+                AB(i, j) += A(i, k)*Detail::conj(BT(j, k));
             }
         }
     }

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -53,6 +53,7 @@ SourceFiles
 #include "Pair.H"
 #include "Field.H"
 #include "autoPtr.H"
+#include "complex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,7 +75,7 @@ class Matrix
 {
     // Private Data
 
-        //- Number of rows and columns in Matrix.
+        //- Number of rows and columns in Matrix
         label mRows_, nCols_;
 
         //- Vector of values of type Type
@@ -86,11 +87,11 @@ class Matrix
         //- Allocate storage for the contents
         inline void doAlloc();
 
-        //- Multiply matrix with vector (A * x)
+        //- Right-multiply Matrix by a column vector (A * x)
         template<class ListType>
         tmp<Field<Type>> AmulImpl(const ListType& x) const;
 
-        //- Multiply matrix transpose with vector (AT * x, or x * A)
+        //- Left-multiply Matrix by a row vector (x * A)
         template<class ListType>
         tmp<Field<Type>> TmulImpl(const ListType& x) const;
 
@@ -138,7 +139,7 @@ public:
         //- Construct given number of rows/columns
         inline explicit Matrix(const labelPair& dims);
 
-        //- Construct given number of rows/columns.
+        //- Construct given number of rows/columns
         //- initializing all elements to zero
         inline Matrix(const labelPair& dims, const zero);
 
@@ -152,15 +153,15 @@ public:
         //- Move construct
         Matrix(Matrix<Form, Type>&& mat);
 
-        //- Copy constructor from matrix of a different form
+        //- Copy constructor from Matrix of a different form
         template<class Form2>
         explicit Matrix(const Matrix<Form2, Type>& mat);
 
-        //- Construct from a block of another matrix
+        //- Construct from a block of another Matrix
         template<class MatrixType>
         Matrix(const ConstMatrixBlock<MatrixType>& Mb);
 
-        //- Construct from a block of another matrix
+        //- Construct from a block of another Matrix
         template<class MatrixType>
         Matrix(const MatrixBlock<MatrixType>& Mb);
 
@@ -179,34 +180,34 @@ public:
 
     // Access
 
-        //- Return the number of rows.
+        //- Return the number of rows
         inline label m() const noexcept;
 
-        //- Return the number of columns.
+        //- Return the number of columns
         inline label n() const noexcept;
 
-        //- Return the number of elements in matrix (m*n)
+        //- 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)
+        //- Return true if Matrix is empty (i.e., size() is zero)
         inline bool empty() const noexcept;
 
         //- Return const pointer to the first data element, which can also
-        //- be used to address into the matrix contents
+        //- be used to address into Matrix contents
         inline const Type* cdata() const noexcept;
 
         //- Return pointer to the first data element, which can also
-        //- be used to address into the matrix contents
+        //- be used to address into Matrix contents
         inline Type* data() noexcept;
 
-        //- Return const pointer to data in the specified row.
+        //- Return const pointer to data in the specified row
         //  Subscript checking only with FULLDEBUG
         inline const Type* rowData(const label irow) const;
 
-        //- Return pointer to data in the specified row.
+        //- Return pointer to data in the specified row
         //  Subscript checking only with FULLDEBUG
         inline Type* rowData(const label irow);
 
@@ -219,87 +220,108 @@ public:
         inline Type& at(const label idx);
 
 
+    // Block Access (const)
 
-    // Block Access
+        //- Return const column or column's subset of Matrix
-
+        //  Return entire column by its index: M.subColumn(a);
-            inline ConstMatrixBlock<mType> block
+        //  Return subset of a column starting from rowIndex: M.subColumn(a, b);
+        //  Return subset of a column starting from rowIndex with szRows elems:
+        //      M.subColumn(a, b, c);
+        inline ConstMatrixBlock<mType> subColumn
         (
-                const label m,
+            const label colIndex,
-                const label n,
+            const label rowIndex = 0,
-                const label mStart,
+            label len = -1
-                const label nStart
         ) const;
 
+        //- Return const row or const row's subset of Matrix
+        //  Return entire row by its index: M.subRow(a);
+        //  Return subset of a row starting from columnIndex: M.subRow(a,b);
+        //  Return subset of a row starting from columnIndex with szCols elems:
+        //      M.subRow(a, b, c);
+        inline ConstMatrixBlock<mType> subRow
+        (
+            const label rowIndex,
+            const label colIndex = 0,
+            label len = -1
+        ) const;
+
+        //- Return const sub-block of Matrix
+        //  Sub-block starting at columnIndex & rowIndex indices
+        inline ConstMatrixBlock<mType> subMatrix
+        (
+            const label rowIndex,
+            const label colIndex,
+            label szRows = -1,
+            label szCols = -1
+        ) const;
+
+        //- Access Field as a ConstMatrixBlock
         template<class VectorSpace>
         inline ConstMatrixBlock<mType> block
         (
-                const label mStart,
+            const label rowIndex,
-                const label nStart
+            const label colIndex
-            ) const;
-
-            inline ConstMatrixBlock<mType> col
-            (
-                const label m,
-                const label rowStart
-            ) const;
-
-            inline ConstMatrixBlock<mType> col
-            (
-                const label m,
-                const label mStart,
-                const label nStart
         ) const;
 
 
-            inline MatrixBlock<mType> block
+    // Block Access (non-const)
+
+        //- Return column or column's subset of Matrix
+        inline MatrixBlock<mType> subColumn
         (
-                const label m,
+            const label colIndex,
-                const label n,
+            const label rowIndex = 0,
-                const label mStart,
+            label len = -1
-                const label nStart
         );
 
+        //- Return row or row's subset of Matrix
+        inline MatrixBlock<mType> subRow
+        (
+            const label rowIndex,
+            const label colIndex = 0,
+            label len = -1
+        );
+
+        //- Return sub-block of Matrix
+        inline MatrixBlock<mType> subMatrix
+        (
+            const label rowIndex,
+            const label colIndex,
+            label szRows = -1,
+            label szCols = -1
+        );
+
+        //- Access Field as a MatrixBlock
         template<class VectorSpace>
         inline MatrixBlock<mType> block
         (
-                const label mStart,
+            const label rowIndex,
-                const label nStart
+            const label colIndex
         );
 
-            inline MatrixBlock<mType> col
-            (
-                const label m,
-                const label rowStart
-            );
-
-            inline MatrixBlock<mType> col
-            (
-                const label m,
-                const label mStart,
-                const label nStart
-            );
 
     // Check
 
-        //- Check index i is within valid range (0 ... m-1).
+        //- Check index i is within valid range [0, m)
         inline void checki(const label irow) const;
 
-        //- Check index j is within valid range (0 ... n-1).
+        //- Check index j is within valid range [0, n)
         inline void checkj(const label jcol) const;
 
         //- Check that dimensions are positive, non-zero
         inline void checkSize() const;
 
-        //- True if all entries have identical values, and matrix is non-empty
+        //- True if all entries have identical values, and Matrix is non-empty
         inline bool uniform() const;
 
 
     // Edit
 
-        //- Clear the Matrix, i.e. set sizes to zero.
+        //- Clear Matrix, i.e. set sizes to zero
         void clear();
 
-        //- Release storage management of the Matrix contents by transferring
+        //- Release storage management of Matrix contents by transferring
         //- management to a List
         List<Type> release();
 
@@ -307,44 +329,65 @@ public:
         void swap(Matrix<Form, Type>& mat);
 
         //- Transfer the contents of the argument Matrix into this Matrix
-        //- and annul the argument Matrix.
+        //- and annul the argument Matrix
         void transfer(Matrix<Form, Type>& mat);
 
-        //- Change the matrix dimensions, preserving the elements
+        //- Change Matrix dimensions, preserving the elements
         void resize(const label m, const label n);
 
-        //- Change the matrix dimensions, preserving the elements
+        //- Change Matrix dimensions, preserving the elements
         inline void setSize(const label m, const label n);
 
-        //- Resize the matrix without reallocating storage (unsafe)
+        //- Resize Matrix without reallocating storage (unsafe)
         inline void shallowResize(const label m, const label n);
 
+        //- Round to zero elements with magnitude smaller than tol (SMALL)
+        void round(const scalar tol = SMALL);
+
 
     // Operations
 
-        //- Return the transpose of the matrix
+        //- Return (conjugate) transpose of Matrix
         Form T() const;
 
-        //- Multiply matrix with vector (A * x)
+        //- Right-multiply Matrix by a column vector (A * x)
         inline tmp<Field<Type>> Amul(const UList<Type>& x) const;
 
-        //- Multiply matrix with vector (A * x)
+        //- Right-multiply Matrix by a column vector (A * x)
         template<class Addr>
         inline tmp<Field<Type>> Amul
         (
             const IndirectListBase<Type, Addr>& x
         ) const;
 
-        //- Multiply matrix transpose with vector (AT * x, or x * A)
+        //- Left-multiply Matrix by a row vector (x * A)
         inline tmp<Field<Type>> Tmul(const UList<Type>& x) const;
 
-        //- Multiply matrix transpose with vector (AT * x, or x * A)
+        //- Left-multiply Matrix by a row vector (x * A)
         template<class Addr>
         inline tmp<Field<Type>> Tmul
         (
             const IndirectListBase<Type, Addr>& x
         ) const;
 
+        //- Extract the diagonal elements. Method may change in the future.
+        List<Type> diag() const;
+
+        //- Assign diagonal of Matrix
+        void diag(const UList<Type>& list);
+
+        //- Return the trace
+        Type trace() const;
+
+        //- Return L2-Norm of chosen column
+        //  Optional without sqrt for parallel usage.
+        scalar columnNorm(const label colIndex, const bool noSqrt=false) const;
+
+        //- Return Frobenius norm of Matrix
+        //  Optional without sqrt for parallel usage.
+        //  https://en.wikipedia.org/wiki/Matrix_norm#Frobenius_norm
+        scalar norm(const bool noSqrt=false) const;
+
 
     // Member Operators
 
@@ -370,11 +413,11 @@ public:
         //- Move assignment
         void operator=(Matrix<Form, Type>&& mat);
 
-        //- Assignment to a block of another matrix
+        //- Assignment to a block of another Matrix
         template<class MatrixType>
         void operator=(const ConstMatrixBlock<MatrixType>& Mb);
 
-        //- Assignment to a block of another matrix
+        //- Assignment to a block of another Matrix
         template<class MatrixType>
         void operator=(const MatrixBlock<MatrixType>& Mb);
 
@@ -390,14 +433,20 @@ public:
         //- Matrix subtraction
         void operator-=(const Matrix<Form, Type>& other);
 
+        //- Matrix scalar addition
+        void operator+=(const Type& s);
+
+        //- Matrix scalar subtraction
+        void operator-=(const Type& s);
+
         //- Matrix scalar multiplication
-        void operator*=(const scalar s);
+        void operator*=(const Type& s);
 
         //- Matrix scalar division
-        void operator/=(const scalar s);
+        void operator/=(const Type& s);
 
 
-    // Random access iterator (non-const)
+    // Random Access Iterator (non-const)
 
         //- Return an iterator to begin traversing a Matrix
         inline iterator begin();
@@ -406,7 +455,7 @@ public:
         inline iterator end();
 
 
-    // Random access iterator (const)
+    // Random Access Iterator (const)
 
         //- Return const_iterator to begin traversing a constant Matrix
         inline const_iterator cbegin() const;
@@ -447,10 +496,73 @@ public:
         {
             return v_;
         }
+
+        //- Deprecated(2019-04) - use subMatrix()
+        //  \deprecated(2019-04) - use subMatrix()
+        ConstMatrixBlock<mType>
+        FOAM_DEPRECATED_FOR(2019-04, "subMatrix() method") block
+        (
+            const label m,
+            const label n,
+            const label mStart,
+            const label nStart
+        ) const
+        {
+            return this->subMatrix(mStart, nStart, m, n);
+        }
+
+        //- Deprecated(2019-04) - use subMatrix()
+        //  \deprecated(2019-04) - use subMatrix()
+        MatrixBlock<mType>
+        FOAM_DEPRECATED_FOR(2019-04, "subMatrix() method") block
+        (
+            const label m,
+            const label n,
+            const label mStart,
+            const label nStart
+        )
+        {
+            return this->subMatrix(mStart, nStart, m, n);
+        }
+
+
+        //- Deprecated(2019-04) - use subColumn()
+        //  \deprecated(2019-04) - use subColumn()
+        ConstMatrixBlock<mType>
+        FOAM_DEPRECATED_FOR(2019-04, "subColumn() method") col
+        (
+            const label m,
+            const label mStart,
+            const label nStart
+        ) const
+        {
+            return this->subColumn(nStart, mStart, m);
+        }
+
+        //- Deprecated(2019-04) - use subColumn()
+        //  \deprecated(2019-04) - use subColumn()
+        MatrixBlock<mType>
+        FOAM_DEPRECATED_FOR(2019-04, "subColumn() method") col
+        (
+            const label m,
+            const label mStart,
+            const label nStart
+        )
+        {
+            return this->subColumn(nStart, mStart, m);
+        }
+
+        //- Deleted(2019-04) - use subColumn()
+        //  \deprecated(2019-04) - use subColumn()
+        void col(const label m, const label rowStart) const = delete;
+
+        //- Deleted(2019-04) - use subColumn()
+        //  \deprecated(2019-04) - use subColumn()
+        void col(const label m, const label rowStart) = delete;
 };
 
 
-// IOstream Operators
+// * * * * * * * * * * * * * * * IOstream Operator  * * * * * * * * * * * * * //
 
 //- Read Matrix from Istream, discarding contents of existing Matrix.
 template<class Form, class Type>
@@ -462,70 +574,100 @@ template<class Form, class Type>
 Ostream& operator<<(Ostream& os, const Matrix<Form, Type>& mat);
 
 
-// Global Functions, Operators
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
 
-//- Find max value in the matrix
+//- Find max value in Matrix
 template<class Form, class Type>
 const Type& max(const Matrix<Form, Type>& mat);
 
-//- Find min value in the matrix
+//- Find min value in Matrix
 template<class Form, class Type>
 const Type& min(const Matrix<Form, Type>& mat);
 
-//- Find the min/max values of the matrix
+//- Find the min/max values of Matrix
 template<class Form, class Type>
 MinMax<Type> minMax(const Matrix<Form, Type>& mat);
 
+
+// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * * //
+
 //- Matrix negation
 template<class Form, class Type>
 Form operator-(const Matrix<Form, Type>& mat);
 
-//- Matrix addition
+//- Matrix addition. Returns Matrix of the same form as the first parameter.
-template<class Form, class Type>
+template<class Form1, class Form2, class Type>
-Form operator+
+Form1 operator+
 (
-    const Matrix<Form, Type>& A,
+    const Matrix<Form1, Type>& A,
-    const Matrix<Form, Type>& B
+    const Matrix<Form2, Type>& B
 );
 
-
+//- Matrix subtraction. Returns Matrix of the same form as the first parameter.
-//- Matrix subtraction
+template<class Form1, class Form2, class Type>
-template<class Form, class Type>
+Form1 operator-
-Form operator-
 (
-    const Matrix<Form, Type>& A,
+    const Matrix<Form1, Type>& A,
-    const Matrix<Form, Type>& B
+    const Matrix<Form2, Type>& B
 );
 
-
+//- Scalar multiplication of Matrix
-//- Scalar multiplication of a matrix
 template<class Form, class Type>
 Form operator*
 (
-    const scalar s,
+    const Type& s,
     const Matrix<Form, Type>& mat
 );
 
+//- Scalar addition of Matrix
+template<class Form, class Type>
+Form operator+
+(
+    const Type& s,
+    const Matrix<Form, Type>& mat
+);
 
-//- Scalar multiplication of a matrix
+//- Scalar subtraction of Matrix
+template<class Form, class Type>
+Form operator-
+(
+    const Type& s,
+    const Matrix<Form, Type>& mat
+);
+
+//- Scalar multiplication of Matrix
 template<class Form, class Type>
 Form operator*
 (
     const Matrix<Form, Type>& mat,
-    const scalar s
+    const Type& s
 );
 
+//- Scalar addition of Matrix
+template<class Form, class Type>
+Form operator+
+(
+    const Matrix<Form, Type>& mat,
+    const Type& s
+);
 
-//- Scalar division of a matrix
+//- Scalar subtraction of Matrix
+template<class Form, class Type>
+Form operator-
+(
+    const Matrix<Form, Type>& mat,
+    const Type& s
+);
+
+//- Scalar division of Matrix
 template<class Form, class Type>
 Form operator/
 (
     const Matrix<Form, Type>& mat,
-    const scalar s
+    const Type& s
 );
 
-
+//- Matrix-Matrix multiplication using ikj-algorithm
-//- Matrix-matrix multiplication
 template<class Form1, class Form2, class Type>
 typename typeOfInnerProduct<Type, Form1, Form2>::type
 operator*
@@ -534,7 +676,6 @@ operator*
     const Matrix<Form2, Type>& B
 );
 
-
 //- Matrix-vector multiplication (A * x), where x is a column vector
 template<class Form, class Type>
 inline tmp<Field<Type>> operator*
@@ -551,8 +692,7 @@ inline tmp<Field<Type>> operator*
     const IndirectListBase<Type, Addr>& x
 );
 
-
+//- Vector-Matrix multiplication (x * A), where x is a row vector
-//- Vector-matrix multiplication (x * A), where x is a row vector
 template<class Form, class Type>
 inline tmp<Field<Type>> operator*
 (
@@ -560,7 +700,7 @@ inline tmp<Field<Type>> operator*
     const Matrix<Form, Type>& mat
 );
 
-//- Vector-matrix multiplication (x * A), where x is a row vector
+//- Vector-Matrix multiplication (x * A), where x is a row vector
 template<class Form, class Type, class Addr>
 inline tmp<Field<Type>> operator*
 (
@@ -568,6 +708,24 @@ inline tmp<Field<Type>> operator*
     const Matrix<Form, Type>& mat
 );
 
+//- Implicit inner product of Matrix-Matrix, equivalent to A.T()*B
+template<class Form1, class Form2, class Type>
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+operator&
+(
+    const Matrix<Form1, Type>& ATranspose,
+    const Matrix<Form2, Type>& B
+);
+
+//- Implicit outer product of Matrix-Matrix, equivalent to A*B.T()
+template<class Form1, class Form2, class Type>
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+operator^
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& BTranspose
+);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/matrices/Matrix/MatrixI.H

@@ -128,40 +128,36 @@ inline bool Foam::Matrix<Form, Type>::empty() const noexcept
 template<class Form, class Type>
 inline void Foam::Matrix<Form, Type>::checki(const label i) const
 {
-    #ifdef FULLDEBUG
     if (!mRows_ || !nCols_)
     {
         FatalErrorInFunction
             << "Attempt to access element from empty matrix"
             << abort(FatalError);
     }
-    if (i < 0 || i >= mRows_)
+    if (i < 0 || mRows_ <= i)
     {
         FatalErrorInFunction
             << "Index " << i << " out of range 0 ... " << mRows_-1
             << abort(FatalError);
     }
-    #endif
 }
 
 
 template<class Form, class Type>
 inline void Foam::Matrix<Form, Type>::checkj(const label j) const
 {
-    #ifdef FULLDEBUG
     if (!mRows_ || !nCols_)
     {
         FatalErrorInFunction
             << "Attempt to access element from empty matrix"
             << abort(FatalError);
     }
-    if (j < 0 || j >= nCols_)
+    if (j < 0 || nCols_ <= j)
     {
         FatalErrorInFunction
             << "index " << j << " out of range 0 ... " << nCols_-1
             << abort(FatalError);
     }
-    #endif
 }
 
 
@@ -238,7 +234,7 @@ template<class Form, class Type>
 inline const Type& Foam::Matrix<Form, Type>::at(const label idx) const
 {
     #ifdef FULLDEBUG
-    if (idx < 0 || idx >= this->size())
+    if (idx < 0 || this->size() <= idx)
     {
         FatalErrorInFunction
             << "index " << idx << " out of range 0 ... " << this->size()
@@ -253,7 +249,7 @@ template<class Form, class Type>
 inline Type& Foam::Matrix<Form, Type>::at(const label idx)
 {
     #ifdef FULLDEBUG
-    if (idx < 0 || idx >= this->size())
+    if (idx < 0 || this->size() <= idx)
     {
         FatalErrorInFunction
             << "index " << idx << " out of range 0 ... " << this->size()
@@ -266,21 +262,74 @@ inline Type& Foam::Matrix<Form, Type>::at(const label idx)
 
 template<class Form, class Type>
 inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block
+Foam::Matrix<Form, Type>::subColumn
 (
-    const label m,
+    const label colIndex,
-    const label n,
+    const label rowIndex,
-    const label mStart,
+    label len
-    const label nStart
 ) const
 {
+    if (len < 0)
+    {
+        len = mRows_ - rowIndex;
+    }
+
     return ConstMatrixBlock<mType>
     (
         *this,
-        m,
+        len, // rows
-        n,
+        1,
-        mStart,
+        rowIndex,
-        nStart
+        colIndex
+    );
+}
+
+
+template<class Form, class Type>
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
+Foam::Matrix<Form, Type>::subRow
+(
+    const label rowIndex,
+    const label colIndex,
+    label len
+) const
+{
+    if (len < 0)
+    {
+        len = nCols_ - colIndex;
+    }
+
+    return ConstMatrixBlock<mType>
+    (
+        *this,
+        1,
+        len, // columns
+        rowIndex,
+        colIndex
+    );
+}
+
+
+template<class Form, class Type>
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
+Foam::Matrix<Form, Type>::subMatrix
+(
+    const label rowIndex,
+    const label colIndex,
+    label szRows,
+    label szCols
+) const
+{
+    if (szRows < 0) szRows = mRows_ - rowIndex;
+    if (szCols < 0) szCols = nCols_ - colIndex;
+
+    return ConstMatrixBlock<mType>
+    (
+        *this,
+        szRows,
+        szCols,
+        rowIndex,
+        colIndex
     );
 }
 
@@ -290,8 +339,8 @@ template<class VectorSpace>
 inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::block
 (
-    const label mStart,
+    const label rowIndex,
-    const label nStart
+    const label colIndex
 ) const
 {
     return ConstMatrixBlock<mType>
@@ -299,68 +348,82 @@ Foam::Matrix<Form, Type>::block
         *this,
         VectorSpace::mRows,
         VectorSpace::nCols,
-        mStart,
+        rowIndex,
-        nStart
+        colIndex
-    );
-}
-
-
-template<class Form, class Type>
-inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col
-(
-    const label m,
-    const label mStart
-) const
-{
-    return ConstMatrixBlock<mType>
-    (
-        *this,
-        m,
-        1,
-        mStart,
-        0
-    );
-}
-
-
-template<class Form, class Type>
-inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col
-(
-    const label m,
-    const label mStart,
-    const label nStart
-) const
-{
-    return ConstMatrixBlock<mType>
-    (
-        *this,
-        m,
-        1,
-        mStart,
-        nStart
     );
 }
 
 
 template<class Form, class Type>
 inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block
+Foam::Matrix<Form, Type>::subColumn
 (
-    const label m,
+    const label colIndex,
-    const label n,
+    const label rowIndex,
-    const label mStart,
+    label len
-    const label nStart
 )
 {
+    if (len < 0)
+    {
+        len = mRows_ - rowIndex;
+    }
+
     return MatrixBlock<mType>
     (
         *this,
-        m,
+        len, // rows
-        n,
+        1,
-        mStart,
+        rowIndex,
-        nStart
+        colIndex
+    );
+}
+
+
+template<class Form, class Type>
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
+Foam::Matrix<Form, Type>::subRow
+(
+    const label rowIndex,
+    const label colIndex,
+    label len
+)
+{
+    if (len < 0)
+    {
+        len = nCols_ - colIndex;
+    }
+
+    return MatrixBlock<mType>
+    (
+        *this,
+        1,
+        len, // columns
+        rowIndex,
+        colIndex
+    );
+}
+
+
+template<class Form, class Type>
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
+Foam::Matrix<Form, Type>::subMatrix
+(
+    const label rowIndex,
+    const label colIndex,
+    label szRows,
+    label szCols
+)
+{
+    if (szRows < 0) szRows = mRows_ - rowIndex;
+    if (szCols < 0) szCols = nCols_ - colIndex;
+
+    return MatrixBlock<mType>
+    (
+        *this,
+        szRows,
+        szCols,
+        rowIndex,
+        colIndex
     );
 }
 
@@ -368,50 +431,19 @@ Foam::Matrix<Form, Type>::block
 template<class Form, class Type>
 template<class VectorSpace>
 inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
+Foam::Matrix<Form, Type>::block
+(
+    const label rowIndex,
+    const label colIndex
+)
 {
     return MatrixBlock<mType>
     (
         *this,
         VectorSpace::mRows,
         VectorSpace::nCols,
-        mStart,
+        rowIndex,
-        nStart
+        colIndex
-    );
-}
-
-
-template<class Form, class Type>
-inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col(const label m, const label mStart)
-{
-    return MatrixBlock<mType>
-    (
-        *this,
-        m,
-        1,
-        mStart,
-        0
-    );
-}
-
-
-template<class Form, class Type>
-inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col
-(
-    const label m,
-    const label mStart,
-    const label nStart
-)
-{
-    return MatrixBlock<mType>
-    (
-        *this,
-        m,
-        1,
-        mStart,
-        nStart
     );
 }
 
@@ -532,8 +564,10 @@ inline const Type& Foam::Matrix<Form, Type>::operator()
     const label jcol
 ) const
 {
+    #ifdef FULLDEBUG
     checki(irow);
     checkj(jcol);
+    #endif
     return v_[irow*nCols_ + jcol];
 }
 
@@ -545,8 +579,10 @@ inline Type& Foam::Matrix<Form, Type>::operator()
     const label jcol
 )
 {
+    #ifdef FULLDEBUG
     checki(irow);
     checkj(jcol);
+    #endif
     return v_[irow*nCols_ + jcol];
 }
 
@@ -554,7 +590,9 @@ inline Type& Foam::Matrix<Form, Type>::operator()
 template<class Form, class Type>
 inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
 {
+    #ifdef FULLDEBUG
     checki(irow);
+    #endif
     return v_ + irow*nCols_;
 }
 
@@ -562,7 +600,9 @@ inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
 template<class Form, class Type>
 inline Type* Foam::Matrix<Form, Type>::operator[](const label irow)
 {
+    #ifdef FULLDEBUG
     checki(irow);
+    #endif
     return v_ + irow*nCols_;
 }
 

src/OpenFOAM/matrices/Matrix/MatrixTools.C

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MatrixTools.H"
+
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
+
+template<class Form1, class Form2, class Type>
+bool Foam::MatrixTools::equal
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& B,
+    const bool verbose,
+    const scalar relTol,
+    const scalar absTol
+)
+{
+    const label len = A.size();
+
+    if (len != B.size())
+    {
+        if (verbose)
+        {
+            Info<< "Matrices have different sizes: "
+                << len << " vs " << B.size() << nl;
+        }
+        return false;
+    }
+
+    auto iter1 = A.cbegin();
+    auto iter2 = B.cbegin();
+
+    for (label i = 0; i < len; ++i)
+    {
+        if ((absTol + relTol*mag(*iter2)) < Foam::mag(*iter1 - *iter2))
+        {
+            if (verbose)
+            {
+                Info<< "Matrix element " << i
+                    << " differs beyond tolerance: "
+                    << *iter1 << " vs " << *iter2 << nl;
+            }
+            return false;
+        }
+
+        ++iter1;
+        ++iter2;
+    }
+
+    if (verbose)
+    {
+        Info<< "All elements equal within the tolerances" << nl;
+    }
+
+    return true;
+}
+
+
+template<class Container>
+Foam::Ostream& Foam::MatrixTools::printMatrix
+(
+    Ostream& os,
+    const Container& mat
+)
+{
+    os  << mat.m() << ' ' << mat.n();
+
+    if (mat.m() == 1)
+    {
+        // row-vector
+        os  << " (";
+        for (label j = 0; j < mat.n(); ++j)
+        {
+            if (j) os  << ' ';
+            os  << mat(0,j);
+        }
+        os  << ')' << nl;
+    }
+    else if (mat.n() == 1)
+    {
+        // col-vector
+
+        os  << " (";
+        for (label i = 0; i < mat.m(); ++i)
+        {
+            if (i) os  << ' ';
+            os  << mat(i,0);
+        }
+        os  << ')' << nl;
+    }
+    else
+    {
+        // Regular
+
+        os  << nl << '(' << nl;
+
+        for (label i = 0; i < mat.m(); ++i)
+        {
+            os  << '(';
+            for (label j = 0; j < mat.n(); ++j)
+            {
+                if (j) os  << ' ';
+                os  << mat(i,j);
+            }
+            os  << ')' << nl;
+        }
+        os  << ')' << nl;
+    }
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/MatrixTools.H

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::MatrixTools
+
+Description
+    Collection of static functions for matrix-related verifications.
+
+SourceFiles
+    MatrixTools.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MatrixTools_H
+#define MatrixTools_H
+
+#include "Matrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class Ostream;
+
+
+/*---------------------------------------------------------------------------*\
+                       Namespace MatrixTools Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace MatrixTools
+{
+
+//- Compare matrix elements for absolute or relative equality
+template<class Form1, class Form2, class Type>
+bool equal
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& B,
+    const bool verbose = false,
+    const scalar relTol = 1e-5,
+    const scalar absTol = 1e-8
+);
+
+
+//- Simple ASCII output of Matrix, MatrixBlock
+template<class Container>
+Ostream& printMatrix(Ostream& os, const Container& mat);
+
+
+} // End namespace MatrixTools
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "MatrixTools.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2016 OpenFOAM Foundation
@@ -71,6 +71,76 @@ Foam::MatrixBlock<MatrixType>::operator Field<cmptType>() const
 }
 
 
+template<class MatrixType>
+Foam::label Foam::ConstMatrixBlock<MatrixType>::disallow
+(
+    const char* what
+) const
+{
+    FatalErrorInFunction
+        << "Block addresses " << what
+        << " outside matrix or invalid matrix components"
+        << abort(FatalError);
+    return 0;
+}
+
+
+template<class MatrixType>
+Foam::label Foam::MatrixBlock<MatrixType>::disallow
+(
+    const char* what
+) const
+{
+    FatalErrorInFunction
+        << "Block addresses " << what
+        << " outside matrix or invalid matrix components"
+        << abort(FatalError);
+    return 0;
+}
+
+
+template<class MatrixType> void Foam::ConstMatrixBlock<MatrixType>::checkIndex
+(
+    const label i,
+    const label j
+) const
+{
+    if (i < 0 || i >= mRows_)
+    {
+        FatalErrorInFunction
+            << "Index " << i << " is out of range 0 ... " << mRows_ - 1
+            << abort(FatalError);
+    }
+    else if (j < 0 || j >= nCols_)
+    {
+        FatalErrorInFunction
+            << "Index " << j << " is out of range 0 ... " << nCols_ - 1
+            << abort(FatalError);
+    }
+}
+
+
+template<class MatrixType> void Foam::MatrixBlock<MatrixType>::checkIndex
+(
+    const label i,
+    const label j
+) const
+{
+    if (i < 0 || i >= mRows_)
+    {
+        FatalErrorInFunction
+            << "Index " << i << " is out of range 0 ... " << mRows_ - 1
+            << abort(FatalError);
+    }
+    else if (j < 0 || j >= nCols_)
+    {
+        FatalErrorInFunction
+            << "Index " << j << " is out of range 0 ... " << nCols_ - 1
+            << abort(FatalError);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class MatrixType>

src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2016 OpenFOAM Foundation
@@ -65,14 +65,21 @@ class ConstMatrixBlock
         //- Const reference to the parent matrix
         const MatrixType& matrix_;
 
-        // Block size
+        //- Block size
         const label mRows_;
         const label nCols_;
 
-        // Block location in parent matrix
+        //- Block location in parent matrix
         const label rowStart_;
         const label colStart_;
 
+
+    // Private Member Functions
+
+        //- Error message for failed sanity checks during matrix construction
+        label disallow(const char *what) const;
+
+
 public:
 
     typedef typename MatrixType::cmptType cmptType;
@@ -110,6 +117,9 @@ public:
 
         //- Convert a column of a matrix to a Field
         operator Field<cmptType>() const;
+
+        //- Check if (i, j) is within range of row-column limits
+        void checkIndex(const label i, const label j) const;
 };
 
 
@@ -125,14 +135,21 @@ class MatrixBlock
         //- Reference to the parent matrix
         MatrixType& matrix_;
 
-        // Block size
+        //- Block size
         const label mRows_;
         const label nCols_;
 
-        // Block location in parent matrix
+        //- Block location in parent matrix
         const label rowStart_;
         const label colStart_;
 
+
+    // Private Member Functions
+
+        //- Error message for failed sanity checks during matrix construction
+        label disallow(const char *what) const;
+
+
 public:
 
     typedef typename MatrixType::cmptType cmptType;
@@ -174,8 +191,11 @@ public:
         //- Convert a column of a matrix to a Field
         operator Field<cmptType>() const;
 
+        //- Check if (i, j) is within range of row-column limits
+        void checkIndex(const label i, const label j) const;
 
-    // Member operators
+
+    // Member Operators
 
         //- Assignment to a compatible matrix
         template<class Form>

src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2016 OpenFOAM Foundation
@@ -38,24 +38,21 @@ Foam::ConstMatrixBlock<MatrixType>::ConstMatrixBlock
 )
 :
     matrix_(matrix),
-    mRows_(m),
+    mRows_(0 < m ? m : disallow("row dim")),
-    nCols_(n),
+    nCols_(0 < n ? n : disallow("col dim")),
-    rowStart_(mStart),
+    rowStart_
-    colStart_(nStart)
-{
-    #ifdef FULLDEBUG
-    if
     (
-        rowStart_ + mRows_ > matrix.m()
+        0 <= mStart
-     || colStart_ + nCols_ > matrix.n()
+     || mStart + mRows_ <= matrix.m()
+     ?  mStart : disallow("row start")
+    ),
+    colStart_
+    (
+        0 <= nStart
+     || nStart + nCols_ <= matrix.n()
+     ?  nStart : disallow("col start")
     )
-    {
+{}
-        FatalErrorInFunction
-            << "Block addresses outside matrix"
-            << abort(FatalError);
-    }
-    #endif
-}
 
 
 template<class MatrixType>
@@ -69,24 +66,21 @@ Foam::MatrixBlock<MatrixType>::MatrixBlock
 )
 :
     matrix_(matrix),
-    mRows_(m),
+    mRows_(0 < m ? m : disallow("row dim")),
-    nCols_(n),
+    nCols_(0 < n ? n : disallow("col dim")),
-    rowStart_(mStart),
+    rowStart_
-    colStart_(nStart)
-{
-    #ifdef FULLDEBUG
-    if
     (
-        rowStart_ + mRows_ > matrix.m()
+        0 <= mStart
-     || colStart_ + nCols_ > matrix.n()
+     || mStart + mRows_ <= matrix.m()
+     ?  mStart : disallow("row start")
+    ),
+    colStart_
+    (
+        0 <= nStart
+     || nStart + nCols_ <= matrix.n()
+     ?  nStart : disallow("col start")
     )
-    {
+{}
-        FatalErrorInFunction
-            << "Block addresses outside matrix"
-            << abort(FatalError);
-    }
-    #endif
-}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -144,18 +138,7 @@ Foam::ConstMatrixBlock<MatrixType>::operator()
 ) const
 {
     #ifdef FULLDEBUG
-    if (i<0 || i>=mRows_)
+    checkIndex(i, j);
-    {
-        FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << mRows_-1
-            << abort(FatalError);
-    }
-    if (j<0 || j>=nCols_)
-    {
-        FatalErrorInFunction
-            << "Index " << j << " out of range 0 ... " << nCols_-1
-            << abort(FatalError);
-    }
     #endif
 
     return matrix_(i + rowStart_, j + colStart_);
@@ -171,18 +154,7 @@ Foam::MatrixBlock<MatrixType>::operator()
 ) const
 {
     #ifdef FULLDEBUG
-    if (i<0 || i>=mRows_)
+    checkIndex(i, j);
-    {
-        FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << mRows_-1
-            << abort(FatalError);
-    }
-    if (j<0 || j>=nCols_)
-    {
-        FatalErrorInFunction
-            << "Index " << j << " out of range 0 ... " << nCols_-1
-            << abort(FatalError);
-    }
     #endif
 
     return matrix_(i + rowStart_, j + colStart_);
@@ -198,18 +170,7 @@ Foam::MatrixBlock<MatrixType>::operator()
 )
 {
     #ifdef FULLDEBUG
-    if (i<0 || i>=mRows_)
+    checkIndex(i, j);
-    {
-        FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << mRows_-1
-            << abort(FatalError);
-    }
-    if (j<0 || j>=nCols_)
-    {
-        FatalErrorInFunction
-            << "Index " << j << " out of range 0 ... " << nCols_-1
-            << abort(FatalError);
-    }
     #endif
 
     return matrix_(i + rowStart_, j + colStart_);

src/OpenFOAM/matrices/QRMatrix/QRMatrix.C

@@ -282,7 +282,7 @@ Foam::QRMatrix<MatrixType>::inv() const
             x[j] = Q_(i, j);
         }
         solvex(x);
-        inv.block(m, 1, 0, i) = x;
+        inv.subColumn(i) = x;
     }
 
     return inv;

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H

@@ -118,7 +118,7 @@ public:
 };
 
 
-// Global functions and operators
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 template<class Type>
 class typeOfInnerProduct<Type, RectangularMatrix<Type>, RectangularMatrix<Type>>
@@ -128,6 +128,7 @@ public:
     typedef RectangularMatrix<Type> type;
 };
 
+
 template<class Type>
 class typeOfInnerProduct<Type, RectangularMatrix<Type>, SquareMatrix<Type>>
 {
@@ -136,6 +137,7 @@ public:
     typedef RectangularMatrix<Type> type;
 };
 
+
 template<class Type>
 class typeOfInnerProduct<Type, SquareMatrix<Type>, RectangularMatrix<Type>>
 {
@@ -145,6 +147,8 @@ public:
 };
 
 
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
 template<class Type>
 RectangularMatrix<Type> outer(const Field<Type>& f1, const Field<Type>& f2);
 

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H

@@ -190,9 +190,9 @@ inline Foam::RectangularMatrix<Type> outer
 {
     RectangularMatrix<Type> f1f2T(f1.size(), f2.size());
 
-    for (label i=0; i<f1f2T.m(); i++)
+    for (label i = 0; i < f1f2T.m(); ++i)
     {
-        for (label j=0; j<f1f2T.n(); j++)
+        for (label j = 0; j < f1f2T.n(); ++j)
         {
             f1f2T(i, j) = f1[i]*f2[j];
         }

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2013-2016 OpenFOAM Foundation
@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SquareMatrix.H"
+#include "RectangularMatrix.H"
 #include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -39,7 +40,7 @@ Foam::scalar Foam::detDecomposed
 {
     Type diagProduct = pTraits<Type>::one;
 
-    for (label i=0; i<matrix.m(); i++)
+    for (label i = 0; i < matrix.m(); ++i)
     {
         diagProduct *= matrix(i, i);
     }

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H

@@ -149,7 +149,7 @@ public:
 };
 
 
-// Global functions and operators
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 //- Return the LU decomposed SquareMatrix det
 template<class Type>
@@ -171,7 +171,6 @@ public:
     typedef SquareMatrix<Type> type;
 };
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H

@@ -251,9 +251,9 @@ inline Foam::SquareMatrix<Type> symmOuter
 {
     SquareMatrix<Type> f1f2T(f1.size());
 
-    for (label i=0; i<f1f2T.m(); i++)
+    for (label i = 0; i < f1f2T.m(); ++i)
     {
-        for (label j=0; j<f1f2T.n(); j++)
+        for (label j = 0; j < f1f2T.n(); ++j)
         {
             f1f2T(i, j) = f1[i]*f2[j];
         }