diff --git a/applications/test/QRMatrix/Make/files b/applications/test/QRMatrix/Make/files
new file mode 100644
index 0000000000..9825bbfcf3
--- /dev/null
+++ b/applications/test/QRMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-QRMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-QRMatrix
diff --git a/applications/test/QRMatrix/Make/options b/applications/test/QRMatrix/Make/options
new file mode 100644
index 0000000000..4e772fdf9d
--- /dev/null
+++ b/applications/test/QRMatrix/Make/options
@@ -0,0 +1,2 @@
+/* EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude */
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/QRMatrix/Test-QRMatrix.C b/applications/test/QRMatrix/Test-QRMatrix.C
new file mode 100644
index 0000000000..fb22d78bde
--- /dev/null
+++ b/applications/test/QRMatrix/Test-QRMatrix.C
@@ -0,0 +1,1484 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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"
+#include "QRMatrix.H"
+#include "Random.H"
+#include "IOmanip.H"
+
+using namespace Foam;
+using namespace Foam::MatrixTools;
+
+#define equal MatrixTools::equal
+#define RUNALL true
+const bool verbose = true;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void horizontalLine()
+{
+    Info<< "+---------+---------+---------+---------+---------+" << nl;
+}
+
+
+// Create random scalar-type matrix
+template<class MatrixType>
+typename std::enable_if
+<
+    !std::is_same<complex, typename MatrixType::cmptType>::value,
+    MatrixType
+>::type makeRandomMatrix
+(
+    const labelPair& dims,
+    Random& rndGen
+)
+{
+    MatrixType mat(dims);
+
+    std::generate
+    (
+        mat.begin(),
+        mat.end(),
+        [&]{ return rndGen.GaussNormal<typename MatrixType::cmptType>(); }
+    );
+
+    return mat;
+}
+
+
+// Create random complex-type matrix
+template<class MatrixType>
+typename std::enable_if
+<
+    std::is_same<complex, typename MatrixType::cmptType>::value,
+    MatrixType
+>::type makeRandomMatrix
+(
+    const labelPair& dims,
+    Random& rndGen
+)
+{
+    MatrixType mat(dims);
+
+    std::generate
+    (
+        mat.begin(),
+        mat.end(),
+        [&]
+        {
+            return complex
+            (
+                rndGen.GaussNormal<scalar>(),
+                rndGen.GaussNormal<scalar>()
+            );
+        }
+    );
+
+    return mat;
+}
+
+
+// Print OpenFOAM matrix in NumPy format
+template<class MatrixType>
+void InfoNumPyFormat
+(
+    const MatrixType& mat,
+    const word objName
+)
+{
+    Info<< objName << ": " << mat.m() << "x" << mat.n() << nl;
+    for (label m = 0; m < mat.m(); ++m)
+    {
+        Info<< "[";
+        for (label n = 0; n < mat.n(); ++n)
+        {
+            if (n == mat.n() - 1)
+            {
+                Info<< mat(m,n);
+            }
+            else
+            {
+                Info<< mat(m,n) << ",";
+            }
+        }
+        if (m == mat.m() - 1)
+        {
+            Info << "]" << nl;
+        }
+        else
+        {
+            Info << "]," << nl;
+        }
+    }
+}
+
+
+// Returns true if two scalars are equal within a given tolerance, and v.v.
+bool isEqual
+(
+    const scalar a,
+    const scalar b,
+    const bool verbose = false,
+    const scalar relTol = 1e-5,
+    const scalar absTol = 1e-8
+)
+{
+    if ((absTol + relTol*mag(b)) < mag(a - b))
+    {
+        if (verbose)
+        {
+            Info<< "Elements are not close in terms of tolerances:"
+                << nl << a << tab << b << nl;
+        }
+        return false;
+    }
+
+    if (verbose)
+    {
+        Info<< "All elems are the same within the tolerances" << nl;
+    }
+
+    return true;
+}
+
+
+// Prints true if a given matrix is empty, and v.v.
+template<class MatrixType>
+void isEmpty
+(
+    const MatrixType& A,
+    const word objName
+)
+{
+    Info<< "Empty " << objName << " = ";
+    if (A.empty())
+    {
+        Info<< "true" << nl;
+    }
+    else
+    {
+        Info<< "false" << nl;
+    }
+}
+
+
+// Checks if Q matrix is a unitary matrix
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+    const MatrixType& Aorig,
+    const MatrixType& Q
+)
+{
+    InfoNumPyFormat(Q, "Q");
+
+    // mathworld.wolfram.com/UnitaryMatrix.html (Retrieved:18-06-19)
+    Info<< nl << "# Q.T()*Q = Q*Q.T() = I:" << nl;
+    const MatrixType QTQ(Q&Q);
+    const MatrixType QQT(Q^Q);
+    MatrixType IMatrix({Q.m(), Q.n()}, Zero);
+    for (label i = 0; i < IMatrix.n(); ++i)
+    {
+        IMatrix(i,i) = pTraits<typename MatrixType::cmptType>::one;
+    }
+    equal(QTQ, IMatrix, verbose);
+    equal(QQT, IMatrix, verbose);
+    equal(QTQ, QQT, verbose);
+}
+
+
+// Checks if R matrix is an upper-triangular matrix
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+    const MatrixType& R
+)
+{
+    InfoNumPyFormat(R, "R");
+
+    // mathworld.wolfram.com/UpperTriangularMatrix.html (Retrieved:18-06-19)
+    Info<< nl << "# R(i, j) = 0 for i > j:" << nl;
+    for (label i = 0; i < R.m(); ++i)
+    {
+        for (label j = 0; j < R.n(); ++j)
+        {
+            if (j < i)
+            {
+                isEqual(mag(R(i, j)), 0.0, verbose);
+            }
+        }
+    }
+}
+
+
+// Checks if the given matrix can be reconstructed by Q and R matrices
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+    const MatrixType& Aorig,
+    const MatrixType& Q,
+    const MatrixType& R
+)
+{
+    InfoNumPyFormat(Aorig, "Aorig");
+
+    cross_check_QRMatrix(Aorig, Q);
+
+    cross_check_QRMatrix(R);
+
+    // mathworld.wolfram.com/QRDecomposition.html (Retrieved:18-06-19)
+    Info<< nl << "# Q*R = A:" << nl;
+    const MatrixType AReconstruct(Q*R);
+    equal(Aorig, AReconstruct, verbose);
+}
+
+
+// Checks if R matrix is an upper-triangular matrix, and obeys column pivoting
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+    const MatrixType& Aorig,
+    const labelList& orderP,
+    const SquareMatrix<typename MatrixType::cmptType>& P,
+    const MatrixType& R
+)
+{
+    InfoNumPyFormat(Aorig, "Aorig");
+
+    InfoNumPyFormat(P, "P");
+
+    cross_check_QRMatrix(R);
+
+    // w.wiki/54m (Retrieved:20-06-19)
+    Info<< nl << "# Diag elems of R non-increasing:"
+        << "|R11| >= |R22| >= ... >= |Rnn|:" << nl;
+    const List<scalar> diag0(mag(R.diag()));
+    List<scalar> diag1(diag0);
+    Foam::sort(diag1, std::greater<scalar>());
+    forAll(diag0, i)
+    {
+        isEqual(diag0[i], diag1[i], verbose);
+    }
+}
+
+
+// Checks if the given matrix can be reconstructed by column-pivoted Q and R
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+    const MatrixType& Aorig,
+    const MatrixType& Q,
+    const labelList& orderP,
+    const SquareMatrix<typename MatrixType::cmptType>& P,
+    const MatrixType& R
+)
+{
+    cross_check_QRMatrix(Aorig, orderP, P, R);
+    cross_check_QRMatrix(Aorig, Q);
+
+    // w.wiki/54m (Retrieved:20-06-19)
+    Info<< nl << "# Q*R = A*P:" << nl;
+    const MatrixType AReconstruct(Q*R);
+    const MatrixType AP(Aorig*P);
+    equal(AP, AReconstruct, verbose);
+}
+
+
+// Checks each constructor of QRMatrix type
+template<class MatrixType>
+void verification_QRMatrix
+(
+    MatrixType& A
+)
+{
+    typedef SquareMatrix<typename MatrixType::cmptType> SMatrix;
+
+    // Create working copies of matrix A
+    const MatrixType Aorig(A);
+
+    // QRMatrix Constructors
+    #if (0 | RUNALL)
+    {
+        QRMatrix<MatrixType> QRNull();
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_R | OUT_OF_PLACE" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+        );
+        QRM.decompose(A);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_R | IN_PLACE" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE
+        );
+        MatrixType A0(A);
+        QRM.decompose(A0);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_QR | OUT_OF_PLACE" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+        );
+        QRM.decompose(A);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        cross_check_QRMatrix(Aorig, Q, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_QR | IN_PLACE" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE
+        );
+        MatrixType A0(A);
+        QRM.decompose(A0);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, Q, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_R | OUT_OF_PLACE | colPivot = on" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        QRM.decompose(A);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(Aorig, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_R | IN_PLACE | colPivot = on" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        MatrixType A0(A);
+        QRM.decompose(A0);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, PList, P, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_QR | OUT_OF_PLACE | colPivot = on" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        QRM.decompose(A);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        cross_check_QRMatrix(Aorig, Q, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# FULL_QR | IN_PLACE | colPivot = on" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        MatrixType A0(A);
+        QRM.decompose(A0);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, Q, PList, P, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_R | OUT_OF_PLACE | colPivot = off" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_R | IN_PLACE | colPivot = off" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_QR | OUT_OF_PLACE | colPivot = off" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        cross_check_QRMatrix(Aorig, Q, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_QR | IN_PLACE | colPivot = off" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, Q, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_R | OUT_OF_PLACE | colPivot = on" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(Aorig, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_R | IN_PLACE | colPivot = on" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, PList, P, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_QR | OUT_OF_PLACE | colPivot = on" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        cross_check_QRMatrix(Aorig, Q, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | FULL_QR | IN_PLACE | colPivot = on" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, Q, PList, P, A0);
+    }
+    #endif
+
+    // constructors with the const type specifier
+    #if (0 | RUNALL)
+    {
+        Info<< "# const A | FULL_R | colPivot = off" << nl;
+        const MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# const A | FULL_QR | colPivot = off" << nl;
+        const MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        cross_check_QRMatrix(Aorig, Q, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# const A | FULL_R | colPivot = on" << nl;
+        const MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_R,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(Aorig, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# const A | FULL_QR | colPivot = on" << nl;
+        const MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        cross_check_QRMatrix(Aorig, Q, PList, P, R);
+    }
+    #endif
+
+    //
+    #if (0 | RUNALL)
+    {
+        Info<< "# REDUCED_R | OUT_OF_PLACE" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+        );
+        QRM.decompose(A);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(R);
+
+        // check if full R and reduced R give the same answer
+        if (A.n() < A.m())
+        {
+            QRMatrix<MatrixType> fullQRM
+            (
+                QRMatrix<MatrixType>::outputTypes::FULL_QR,
+                QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+            );
+            fullQRM.decompose(A);
+            const MatrixType& fullR(fullQRM.R());
+            MatrixType fullR2(fullR.subMatrix(0,0,R.m()));
+            equal(fullR2, R, verbose);
+        }
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# REDUCED_R | IN_PLACE" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE
+        );
+        MatrixType A0(A);
+        QRM.decompose(A0);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(A0);
+
+        // check if full R and reduced R give the same answer
+        if (A.n() < A.m())
+        {
+            QRMatrix<MatrixType> fullQRM
+            (
+                QRMatrix<MatrixType>::outputTypes::FULL_QR,
+                QRMatrix<MatrixType>::storeMethods::IN_PLACE
+            );
+            MatrixType A1(A);
+            fullQRM.decompose(A1);
+            MatrixType fullR(A1.subMatrix(0,0,A0.m()));
+            equal(fullR, A0, verbose);
+        }
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# REDUCED_R | OUT_OF_PLACE | colPivot = on" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        QRM.decompose(A);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(Aorig, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# REDUCED_R | IN_PLACE | colPivot = on" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        MatrixType A0(A);
+        QRM.decompose(A0);
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, PList, P, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | REDUCED_R | OUT_OF_PLACE | colPivot = off" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | REDUCED_R | IN_PLACE | colPivot = off" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | REDUCED_R | OUT_OF_PLACE | colPivot = on" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(Aorig, PList, P, R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# A | REDUCED_R | IN_PLACE | colPivot = on" << nl;
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        isEmpty(R, "R");
+        cross_check_QRMatrix(Aorig, PList, P, A0);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# const A | REDUCED_R | colPivot = off" << nl;
+        const MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::colPivoting::FALSE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(R);
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< "# const A | REDUCED_R | colPivot = on" << nl;
+        const MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+            QRMatrix<MatrixType>::colPivoting::TRUE
+        );
+        const MatrixType& Q(QRM.Q());
+        const MatrixType& R(QRM.R());
+        const labelList& PList(QRM.orderP());
+        const SMatrix P(QRM.P());
+        isEmpty(Q, "Q");
+        cross_check_QRMatrix(Aorig, PList, P, R);
+    }
+    #endif
+}
+
+
+// Checks the direct solution of a linear system by QRMatrix
+template<class MatrixType>
+void verification_QRMatrixSolve
+(
+    MatrixType& A
+)
+{
+    typedef SquareMatrix<typename MatrixType::cmptType> SMatrix;
+
+    // Create working copies of matrix A
+    const MatrixType Aorig(A);
+
+    // Direct solution of a linear system by QR decomposition
+    #if (0 | RUNALL)
+    {
+        MatrixType A0(A);
+        QRMatrix<MatrixType> QRM
+        (
+            A0,
+            QRMatrix<MatrixType>::outputTypes::FULL_QR
+        );
+
+        Info<< nl << "# Solution of A*x = b with A = Q*R:" << nl;
+        const scalarField residual(A0.m(), 0);
+        const scalarField source(A0.m(), 1);
+
+        const scalarField x(QRM.solve(source));
+        const scalarField A0xMinusSource(A0*x - source);
+        const scalarField xMinusQRinvSource(x - QRM.inv()*source);
+        const SMatrix QRinvA0(QRM.inv()*A0);
+        const SMatrix identity(A0.m(), I);
+
+        forAll(residual, i)
+        {
+            isEqual(A0xMinusSource[i], residual[i], verbose);
+            isEqual(xMinusQRinvSource[i], residual[i], verbose);
+        }
+        equal(QRinvA0, identity, verbose);
+    }
+    #endif
+}
+
+
+// Checks the parallel tall-skinny QR decomposition
+template<class Type>
+void verification_tsqr
+(
+    Random& rndGen
+)
+{
+    typedef RectangularMatrix<Type> RMatrix;
+
+    // Size of the full matrix and its partitions
+    const label nColsSum = rndGen.position(1, 100);
+    const label qParts = rndGen.position(10, 30);
+    List<label> mRowsList(qParts, Zero);
+
+    label mRowsSum = 0;
+    for (label i = 0; i < qParts; ++i)
+    {
+        const label mRows = rndGen.position(nColsSum, 10*nColsSum);
+        mRowsList[i] = mRows;
+        mRowsSum += mRows;
+    }
+
+    # if 0
+    Info<< "mRowsSum = " << mRowsSum << nl;
+    Info<< "nColsSum = " << nColsSum << nl;
+    Info<< "mRowsList = " << mRowsList << nl;
+    #endif
+
+    if (mRowsSum < nColsSum)
+    {
+        FatalErrorInFunction
+            << "Tall skinny QR decomposition cannot be executed for matrices "
+            << "with number of columns (nCols) = " << nColsSum << " > "
+            << "than number of rows (mRows) = " << mRowsSum << "."
+            << "Yet nCols > mRows is allowed for a submatrix."
+            << exit(FatalError);
+    }
+
+    // Perform the QR decomposition for the full matrix
+    const RMatrix A(makeRandomMatrix<RMatrix>({mRowsSum, nColsSum}, rndGen));
+    QRMatrix<RMatrix> QRM
+    (
+        A,
+        QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+        QRMatrix<RMatrix>::colPivoting::FALSE
+    );
+    const RMatrix& R = QRM.R();
+
+    RMatrix RMag(R);
+    for (auto& val : RMag)
+    {
+        val = mag(val);
+    }
+
+    // Partition the full matrix,
+    // and perform the QR decomposition on each partition
+    RMatrix subRs({qParts*nColsSum, nColsSum}, Zero);
+    label mRowsIndex = 0;
+    for (label i = 0; i < qParts; ++i)
+    {
+        RMatrix subA({mRowsList[i], nColsSum}, Zero);
+        subA = A.subMatrix(mRowsIndex, 0, mRowsList[i]);
+        mRowsIndex += mRowsList[i];
+
+        QRMatrix<RMatrix> subQRM
+        (
+            subA,
+            QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+            QRMatrix<RMatrix>::colPivoting::FALSE
+        );
+        const RMatrix& subR = subQRM.R();
+
+        // Append subR
+        subRs.subMatrix(i*nColsSum, 0, nColsSum) = subR;
+    }
+
+    // Perform the tall-skinny QR decomposition
+    QRMatrix<RMatrix> TSQR
+    (
+        subRs,
+        QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+        QRMatrix<RMatrix>::colPivoting::FALSE
+    );
+    const RMatrix& TSQRR = TSQR.R();
+
+    RMatrix TSQRRMag(TSQRR);
+    for (auto& val : TSQRRMag)
+    {
+        val = mag(val);
+    }
+
+    // Compare magnitude of R, since R is not unique up to the sign
+    equal(RMag, TSQRRMag, verbose);
+
+    #if 0
+    InfoNumPyFormat(R, "R");
+    InfoNumPyFormat(TSQRR, "TSQRR");
+    #endif
+}
+
+
+// Checks the back substitution function
+template<class Type>
+void verification_backSubstitution
+(
+    const SquareMatrix<Type>& A,
+    const RectangularMatrix<Type>& b
+)
+{
+    // Ax = b
+    typedef RectangularMatrix<Type> RMatrix;
+
+    QRMatrix<RMatrix> QRNull;
+    const RMatrix X(QRNull.backSubstitution(A, b));
+
+    #if 1
+    {
+        InfoNumPyFormat(A, "A");
+        InfoNumPyFormat(b, "b");
+        InfoNumPyFormat(X, "x");
+    }
+    #endif
+
+    #if (0 | RUNALL)
+    {
+        Info<< nl << "# A*X = b:" << nl;
+        const RMatrix AX(A*X);
+        equal(AX, b, verbose, 1e-3, 1e-3);
+    }
+    #endif
+}
+
+
+// Checks the Moore-Penrose pseudo-inverse function
+template<class MatrixType>
+void verification_pinv
+(
+    const MatrixType& A
+)
+{
+    Info<< "### Moore-Penrose pseudo-inverse: pinv()" << nl << nl;
+    const MatrixType APlus(pinv(A));
+
+    #if 0
+    InfoNumPyFormat(A, "A");
+    InfoNumPyFormat(APlus, "A^+");
+    #endif
+
+    if (0 | RUNALL)
+    {
+        Info<< nl << "# (A^+)^+ = A:" << nl;
+        const MatrixType APlusPlus(pinv(APlus));
+        equal(APlusPlus, A, verbose);
+    }
+
+    // The four Moore-Penrose conditions
+    // w.wiki/5RL (Retrieved: 29-06-19)
+    if (0 | RUNALL)
+    {
+        Info<< nl << "# A*(A^+)*A = A:" << nl;
+        const MatrixType AAPlusA((A*APlus)*A);
+        equal(AAPlusA, A, verbose);
+    }
+
+    if (0 | RUNALL)
+    {
+        Info<< nl << "# (A^+)*A*(A^+) = A:" << nl;
+        const MatrixType APlusAAPlus((APlus*A)*APlus);
+        equal(APlusAAPlus, APlus, verbose);
+    }
+
+    if (0 | RUNALL)
+    {
+        Info<< nl << "# (A*(A^+)).T() = A*(A^+):" << nl;
+        const MatrixType AAPlus(A*APlus);
+        const MatrixType AAPlusT(AAPlus.T());
+        equal(AAPlusT, AAPlus, verbose);
+    }
+
+    if (0 | RUNALL)
+    {
+        Info<< nl << "# ((A^+)*A = ((A^+)*A).T():" << nl;
+        const MatrixType APlusA(APlus*A);
+        const MatrixType APlusAT(APlusA.T());
+        equal(APlusA, APlusAT, verbose);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    typedef SquareMatrix<scalar> SMatrix;
+    typedef RectangularMatrix<scalar> RMatrix;
+    typedef SquareMatrix<complex> SCMatrix;
+    typedef RectangularMatrix<complex> RCMatrix;
+
+    Info<< setprecision(15);
+    Random rndGen(1234);
+    label numberOfTests = 10;
+
+    Info<< "### QR Decomposition:" << nl << nl;
+
+    // SquareMatrix<scalar>
+    #if (0 | RUNALL)
+    {
+        horizontalLine();
+
+        Info<< "## " << numberOfTests << " SquareMatrix<scalar> A:" << nl;
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            const label mRows = rndGen.position(1, 50);
+            Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+
+            SMatrix A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
+            #if (0 | RUNALL)
+            {
+                const SMatrix B(A);
+                verification_pinv(B);
+            }
+            #endif
+            verification_QRMatrix(A);
+            verification_QRMatrixSolve(A);
+
+        }
+
+        horizontalLine();
+    }
+    #endif
+
+
+    // RectangularMatrix<scalar>
+    #if (0 | RUNALL)
+    {
+        horizontalLine();
+
+        Info<< "## " << numberOfTests << " RectangularMatrix<scalar> A:" << nl;
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            const label mRows = rndGen.position(1, 50);
+            const label nCols = rndGen.position(1, 50);
+            Info<< nl << nl << "# Random matrix A with"
+                << " random mRows = " << mRows
+                << " random nCols = " << nCols << nl;
+
+            RMatrix A(makeRandomMatrix<RMatrix>({mRows, nCols}, rndGen));
+
+            #if (0 | RUNALL)
+            {
+                const RMatrix B(A);
+                verification_pinv(B);
+            }
+            #endif
+
+            verification_QRMatrix(A);
+        }
+
+        horizontalLine();
+    }
+    #endif
+
+
+    // SquareMatrix<complex(scalar, 0.0)>
+    #if (0 | RUNALL)
+    {
+        horizontalLine();
+
+        Info<< "## " << numberOfTests << " SquareMatrix<complex, 0> A:" << nl;
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            const label mRows = rndGen.position(1, 50);
+            Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+
+            SCMatrix A({mRows, mRows}, complex(0, 0));
+
+            for (auto& val : A)
+            {
+                val.Re() = rndGen.GaussNormal<scalar>();
+            }
+
+            verification_QRMatrix(A);
+        }
+
+        horizontalLine();
+    }
+    #endif
+
+
+    // SquareMatrix<complex(scalar, scalar)>
+    #if (0 | RUNALL)
+    {
+        horizontalLine();
+
+        Info<< "## " << numberOfTests << " SquareMatrix<complex> A:" << nl;
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            const label mRows = rndGen.position(1, 50);
+            Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+
+            SCMatrix A(makeRandomMatrix<SCMatrix>({mRows, mRows}, rndGen));
+
+            #if (0 | RUNALL)
+            {
+                const SCMatrix B(A);
+                verification_pinv(B);
+            }
+            #endif
+
+            verification_QRMatrix(A);
+        }
+
+        horizontalLine();
+    }
+    #endif
+
+
+    // RectangularMatrix<complex(scalar, scalar)>
+    #if (0 | RUNALL)
+    {
+        horizontalLine();
+
+        Info<< "## " << numberOfTests << " RectangularMatrix<complex> A:" << nl;
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            const label mRows = rndGen.position(1, 50);
+            const label nCols = rndGen.position(1, 50);
+            Info<< nl << nl << "# Random matrix A with"
+                << " random mRows = " << mRows
+                << " random nCols = " << nCols << nl;
+
+            RCMatrix A(makeRandomMatrix<RCMatrix>({mRows, nCols}, rndGen));
+
+            #if (0 | RUNALL)
+            {
+                const RCMatrix B(A);
+                verification_pinv(B);
+            }
+            #endif
+
+            verification_QRMatrix(A);
+        }
+
+        horizontalLine();
+    }
+    #endif
+
+
+    #if (0 | RUNALL)
+    {
+        horizontalLine();
+
+        Info<< "### Parallel direct tall-skinny QR decomposition:" << nl;
+        /*
+        Parallel direct tall-skinny QR decomposition:
+            Benson, A. R., Gleich, D. F., & Demmel, J. (2013).
+            Direct QR factorizations for tall-and-skinny matrices in MapReduce
+            architectures.
+            2013 IEEE International Conference on Big Data.
+            DOI:10.1109/bigdata.2013.6691583
+
+            Demmel, J., Grigori, L., Hoemmen, M., & Langou, J. (2012).
+            Communication-optimal parallel and sequential QR and LU
+            factorizations.
+            SIAM Journal on Scientific Computing, 34(1), A206-A239.
+            DOI:10.1137/080731992
+        */
+        // (Benson et al. 2013, Fig. 5) & (Demmel et al. 2012, Fig. 2)
+
+        Info<< "## " << numberOfTests << " <scalar> tests:" << nl;
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            verification_tsqr<scalar>(rndGen);
+        }
+
+        Info<< nl << "## " << numberOfTests << " <complex> tests:" << nl;
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            verification_tsqr<complex>(rndGen);
+        }
+
+        horizontalLine();
+    }
+    #endif
+
+
+    #if (0 | RUNALL)
+    {
+        Info<< "### Back substitution:" << nl << nl;
+
+        numberOfTests = 100;
+
+        // <scalar>
+        #if (0 | RUNALL)
+        {
+            horizontalLine();
+
+            Info<< "## " << numberOfTests << " <scalar>:" << nl;
+
+            for (label i = 0; i < numberOfTests; ++i)
+            {
+                const label mRows = rndGen.position(20, 50);
+                const label pCols = rndGen.position(20, 50);
+                Info<< nl << nl << "# Random square matrix A with"
+                    << " random mRows = " << mRows
+                    << " random nCols = " << mRows << nl;
+
+                SMatrix A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
+
+                // Zeroise the lower diagonal,
+                // so that A becomes an upper-triangular matrix
+                for (label i = 0; i < mRows; ++i)
+                {
+                    for (label j = 0; j < i; ++j)
+                    {
+                        A(i, j) = 0.0;
+                    }
+                }
+
+                // det(triangularA) = prod(diag(triangularA))
+                scalar det = 1;
+                const List<scalar> diag0(A.diag());
+                for (const auto& val : diag0)
+                {
+                    det *= val;
+                }
+
+                if (1e-8 < det)
+                {
+                    Info<< "# Random matrix b with"
+                        << " random mRows = " << mRows
+                        << " random nCols = " << pCols << nl;
+
+                    const RMatrix b
+                    (
+                        makeRandomMatrix<RMatrix>({mRows, pCols}, rndGen)
+                    );
+                    Info<< "det(A) = " << det << nl;
+                    verification_backSubstitution(A, b);
+                }
+                else
+                {
+                    Info<< "Random matrix A is a nearly-singular matrix."
+                        << " Skipping back-substitution" << nl;
+                }
+            }
+
+            horizontalLine();
+        }
+        #endif
+
+        // <complex>
+        #if (0 | RUNALL)
+        {
+            horizontalLine();
+
+            Info<< "## " << numberOfTests << " <complex>:" << nl;
+
+            for (label i = 0; i < numberOfTests; ++i)
+            {
+                const label mRows = rndGen.position(20, 50);
+                const label pCols = rndGen.position(20, 50);
+                Info<< nl << nl << "# Random square matrix A with"
+                    << " random mRows = " << mRows
+                    << " random nCols = " << mRows << nl;
+
+                SCMatrix A(makeRandomMatrix<SCMatrix>({mRows, mRows}, rndGen));
+
+                // Zeroise the lower diagonal,
+                // so that A becomes an upper-triangular matrix
+                for (label i = 0; i < mRows; ++i)
+                {
+                    for (label j = 0; j < i; ++j)
+                    {
+                        A(i, j) = pTraits<complex>::zero;
+                    }
+                }
+
+                // det(triangularA) = prod(diag(triangularA))
+                complex det = pTraits<complex>::one;
+                const List<complex> diag0(A.diag());
+                for (const auto& val : diag0)
+                {
+                    det *= val;
+                }
+
+                if (1e-8 < mag(det))
+                {
+                    Info<< "# Random matrix b with"
+                        << " random mRows = " << mRows
+                        << " random nCols = " << pCols << nl;
+
+                    const RCMatrix b
+                    (
+                        makeRandomMatrix<RCMatrix>({mRows, pCols}, rndGen)
+                    );
+                    Info<< "det(A) = " << det << nl;
+                    verification_backSubstitution(A, b);
+                }
+                else
+                {
+                    Info<< "Random matrix A is a nearly-singular matrix."
+                        << " Skipping back-substitution" << nl;
+                }
+            }
+
+            horizontalLine();
+        }
+        #endif
+    }
+    #endif
+
+
+    Info<< nl << "End" << nl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index f4e5fbb996..4394700e0e 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -28,132 +28,99 @@ License
 
 #include "QRMatrix.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class MatrixType>
-Foam::QRMatrix<MatrixType>::QRMatrix(const MatrixType& mat)
+void Foam::QRMatrix<MatrixType>::qr
+(
+    MatrixType& A
+)
 {
-    decompose(mat);
+    const label nIter = min(A.m() - 1, A.n());
+
+    // Loop through all subcolumns of which the diagonal elem is the first elem
+    for (label k = 0; k < nIter; ++k)
+    {
+        const RMatrix u(householderReflector(A.subColumn(k, k)));
+
+        applyHouseholder(A, u, k);
+    }
+
+    if (outputType_ == outputTypes::REDUCED_R)
+    {
+        A.resize(A.n(), A.n());
+    }
 }
 
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
 template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::decompose(const MatrixType& mat)
+void Foam::QRMatrix<MatrixType>::qrPivot
+(
+    MatrixType& A
+)
 {
-    const label m = mat.m();
-    const label n = mat.n();
+    const label nCols = A.n();
+    const label nIter = min(A.m() - 1, A.n());
 
-    // Initialize the R-matrix to M
-    R_ = mat;
+    // Initialise permutation vector, and column norm vector
+    P_ = identity(nCols);
 
-    // Initialize the Q-matrix to I
-    Q_.setSize(m);
-    Q_ = I;
-
-    // Pre-allocate temporary storage for the Householder steps
-    QMatrixType Qk(m);
-    QMatrixType Rk(m);
-    Field<cmptType> uk(m);
-
-    for (label k=0; k<n; k++)
+    // Initialise vector norms of each column of A
+    List<scalar> colNorms(nCols);
+    for (label k = 0; k < nCols; ++k)
     {
-        // alpha = -|column k of Rk|
-        cmptType alpha = Zero;
-        for (label bi=k; bi<m; bi++)
-        {
-            alpha += sqr(R_(bi, k));
-        }
-        alpha = sqrt(alpha);
+        colNorms[k] = A.columnNorm(k, true);
+    }
 
-        if (R_(k, k) > 0)
+    // Loop through all subcolumns of which the diagonal elem is the first elem
+    for (label k = 0; k < nIter; ++k)
+    {
+        const labelRange colRange(k, nCols);
+        const SubList<scalar> subColNorms(colNorms, colRange);
+
+        // Column pivoting
+        const label maxColNormi =
+            std::distance
+            (
+                subColNorms.cbegin(),
+                std::max_element(subColNorms.cbegin(), subColNorms.cend())
+            );
+
+        // Swap R_, P_ and colNorms_ according to pivot column if the current
+        // column is not the max norm column by avoiding any column swap where
+        // the leading elem is very small
+        if (maxColNormi != 0 && SMALL < mag(A(k, k + maxColNormi)))
         {
-            alpha = -alpha;
+            const RMatrix R1(A.subColumn(k));
+            const RMatrix R2(A.subColumn(maxColNormi + k));
+            A.subColumn(k) = R2;
+            A.subColumn(maxColNormi + k) = R1;
+
+            Swap(P_[k], P_[maxColNormi + k]);
+            Swap(colNorms[k], colNorms[maxColNormi + k]);
         }
 
-        // uk = column k of Rk - alpha*ek
-        // rSumSqrUk = 2/sum(sqr(uk))
-        uk[k] = R_(k, k) - alpha;
-        cmptType rSumSqrUk = sqr(uk[k]);
-        for (label bi=k+1; bi<m; bi++)
         {
-            uk[bi] = R_(bi, k);
-            rSumSqrUk += sqr(uk[bi]);
-        }
-        rSumSqrUk = 2/rSumSqrUk;
+            const RMatrix u(householderReflector(A.subColumn(k, k)));
 
-        // Qk = I - 2*u*uT/sum(sqr(uk))
-        for (label bi=k; bi<m; bi++)
-        {
-            for (label bj=k; bj<m; bj++)
-            {
-                Qk(bi, bj) = -rSumSqrUk*uk[bi]*uk[bj];
-            }
-            Qk(bi, bi) += 1;
+            applyHouseholder(A, u, k);
         }
 
-        // Rk = Qk*R
-        for (label bi=k; bi<m; bi++)
+        // Update norms
+        if (k < nIter - 1)
         {
-            for (label bk=k; bk<n; bk++)
+            label q = k + 1;
+            for (const auto& val : RMatrix(A.subRow(k, q)))
             {
-                Rk(bi, bk) = Zero;
-                for (label bj=k; bj<n; bj++)
-                {
-                    Rk(bi, bk) += Qk(bi, bj)*R_(bj, bk);
-                }
+                colNorms[q] -= magSqr(val);
+                ++q;
             }
         }
+    }
 
-        // Ensure diagonal is positive
-        if (R_(k, k) < 0)
-        {
-            // R = -Rk
-            // Q = -Q
-            for (label bi=k; bi<m; bi++)
-            {
-                for (label bj=k; bj<n; bj++)
-                {
-                    R_(bi, bj) = -Rk(bi, bj);
-                }
-                for (label bj=k; bj<m; bj++)
-                {
-                    Q_(bi, bj) = -Q_(bi, bj);
-                }
-            }
-        }
-        else
-        {
-            // R = Rk
-            for (label bi=k; bi<m; bi++)
-            {
-                for (label bj=k; bj<n; bj++)
-                {
-                    R_(bi, bj) = Rk(bi, bj);
-                }
-            }
-        }
-
-        // Q = Q*Qk (using Rk as temporary storage)
-        for (label bi=0; bi<m; bi++)
-        {
-            for (label bk=k; bk<m; bk++)
-            {
-                Rk(bi, bk) = Zero;
-                for (label bj=k; bj<m; bj++)
-                {
-                    Rk(bi, bk) += Q_(bi, bj)*Qk(bj, bk);
-                }
-            }
-        }
-        for (label bi=0; bi<m; bi++)
-        {
-            for (label bj=k; bj<n; bj++)
-            {
-                Q_(bi, bj) = Rk(bi, bj);
-            }
-        }
+    if (outputType_ == outputTypes::REDUCED_R)
+    {
+        A.resize(A.n(), A.n());
     }
 }
 
@@ -167,7 +134,7 @@ void Foam::QRMatrix<MatrixType>::solvex
 {
     const label n = R_.n();
 
-    for (label i = n - 1; i >= 0; --i)
+    for (label i = n - 1; 0 <= i; --i)
     {
         cmptType sum = x[i];
 
@@ -199,7 +166,7 @@ void Foam::QRMatrix<MatrixType>::solveImpl
     // Assert (&x != &source) ?
     const label m = Q_.m();
 
-    // x = Q_.T()*source;  ie, Q_.Tmul(source)
+    // x = Q_.T()*source;  i.e., Q_.Tmul(source)
     for (label i = 0; i < m; ++i)
     {
         x[i] = 0;
@@ -213,6 +180,179 @@ void Foam::QRMatrix<MatrixType>::solveImpl
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix()
+:
+    outputType_(),
+    storeMethod_(),
+    colPivot_()
+{}
+
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix
+(
+    const outputTypes outputType,
+    const storeMethods storeMethod,
+    const colPivoting colPivot
+)
+:
+    outputType_(outputType),
+    storeMethod_(storeMethod),
+    colPivot_(colPivot),
+    Q_(),
+    R_(),
+    P_()
+{}
+
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix
+(
+    MatrixType& A,
+    const outputTypes outputType,
+    const storeMethods storeMethod,
+    const colPivoting colPivot
+)
+:
+    outputType_(outputType),
+    storeMethod_(storeMethod),
+    colPivot_(colPivot),
+    Q_(),
+    R_(),
+    P_()
+{
+    decompose(A);
+}
+
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix
+(
+    const MatrixType& A,
+    const outputTypes outputType,
+    const colPivoting colPivot
+)
+:
+    outputType_(outputType),
+    storeMethod_(storeMethods::OUT_OF_PLACE),
+    colPivot_(colPivot),
+    Q_(),
+    R_(),
+    P_()
+{
+    decompose(A);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose
+(
+    MatrixType& A
+)
+{
+    // Check whether settings and input are consistent for reduced QRMatrix
+    if (A.m() <= A.n() && outputType_ == outputTypes::REDUCED_R)
+    {
+        outputType_ = outputTypes::FULL_R;
+
+        #if FULLDEBUG
+        WarningInFunction
+            << "Reduced QR decomposition is by definition limited to matrices"
+            << " wherein rows > columns, thus computing FULL decomposition."
+            << nl;
+        #endif
+    }
+
+    // Allocate resources for Q_, if need be
+    if (outputType_ == outputTypes::FULL_QR)
+    {
+        Q_ = MatrixType({A.m(), A.m()}, I);
+    }
+
+    // Allocate resources for R_ and execute decomposition
+    switch (storeMethod_)
+    {
+        case storeMethods::IN_PLACE:
+        {
+            if (colPivot_)
+            {
+                qrPivot(A);
+            }
+            else
+            {
+                qr(A);
+            }
+            break;
+        }
+
+        case storeMethods::OUT_OF_PLACE:
+        {
+            R_ = A;
+
+            if (colPivot_)
+            {
+                qrPivot(R_);
+            }
+            else
+            {
+                qr(R_);
+            }
+            break;
+        }
+    }
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose
+(
+    const MatrixType& A
+)
+{
+    if (storeMethod_ == storeMethods::IN_PLACE)
+    {
+        WarningInFunction
+            << "const type qualifier invalidates storeMethods::IN_PLACE." << nl;
+    }
+
+    // Check whether settings and input are consistent for reduced QRMatrix
+    if (A.m() <= A.n() && outputType_ == outputTypes::REDUCED_R)
+    {
+        outputType_ = outputTypes::FULL_R;
+
+        #if FULLDEBUG
+        WarningInFunction
+            << "Reduced QR decomposition is by definition limited to matrices"
+            << " wherein rows > columns, thus computing FULL decomposition."
+            << nl;
+        #endif
+    }
+
+    // Allocate resources for Q_, if need be
+    if (outputType_ == outputTypes::FULL_QR)
+    {
+        Q_ = MatrixType({A.m(), A.m()}, I);
+    }
+
+    // Allocate resources for R_ and execute decomposition
+    R_ = A;
+
+    if (colPivot_)
+    {
+        qrPivot(R_);
+    }
+    else
+    {
+        qr(R_);
+    }
+}
+
+
 template<class MatrixType>
 void Foam::QRMatrix<MatrixType>::solve
 (
@@ -268,13 +408,13 @@ Foam::QRMatrix<MatrixType>::solve
 
 
 template<class MatrixType>
-typename Foam::QRMatrix<MatrixType>::QMatrixType
+typename Foam::QRMatrix<MatrixType>::SMatrix
 Foam::QRMatrix<MatrixType>::inv() const
 {
     const label m = Q_.m();
 
     Field<cmptType> x(m);
-    QMatrixType inv(m);
+    SMatrix inv(m);
 
     for (label i = 0; i < m; ++i)
     {
@@ -290,4 +430,168 @@ Foam::QRMatrix<MatrixType>::inv() const
 }
 
 
+template<class MatrixType>
+Foam::RectangularMatrix<typename MatrixType::cmptType>
+Foam::QRMatrix<MatrixType>::backSubstitution
+(
+    const SMatrix& A,
+    const RMatrix& rhs
+)
+{
+    const label mRows = A.m();
+    const label nCols = A.n();
+    const label pCols = rhs.n();
+
+    #ifdef FULLDEBUG
+    {
+        const label qRows = rhs.m();
+        if (mRows != qRows)
+        {
+            FatalErrorInFunction
+                << "Linear system is not solvable since the number of rows of"
+                << "A and rhs are not equal:" << tab << mRows << "vs." << qRows
+                << abort(FatalError);
+        }
+
+        const List<cmptType> diag(A.diag());
+        for (const auto& val : diag)
+        {
+            if (mag(val) < SMALL)
+            {
+                WarningInFunction
+                    << "SMALL-valued diagonal elem in back-substitution." << nl;
+            }
+        }
+    }
+    #endif
+
+    RMatrix X({nCols, pCols}, Zero);
+
+    for (label i = 0; i < pCols; ++i)
+    {
+        for (label j = mRows - 1; -1 < j; --j)
+        {
+            cmptType alpha(rhs(j, i));
+
+            for (label k = j + 1; k < mRows; ++k)
+            {
+                alpha -= X(k, i)*A(j, k);
+            }
+
+            X(j, i) = alpha/A(j, j);
+        }
+    }
+
+    return X;
+}
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<class MatrixType>
+MatrixType Foam::pinv
+(
+    const MatrixType& A,
+    const scalar tol
+)
+{
+    typedef typename MatrixType::cmptType cmptType;
+
+    #if FULLDEBUG
+    if (A.empty())
+    {
+        FatalErrorInFunction
+            << "Empty matrix." << abort(FatalError);
+    }
+
+    // Check if A is full-rank
+    #endif
+
+    QRMatrix<MatrixType> QRM
+    (
+        A,
+        QRMatrix<MatrixType>::outputTypes::FULL_QR,
+        QRMatrix<MatrixType>::colPivoting::TRUE
+    );
+    const MatrixType& R(QRM.R());
+    const MatrixType& Q(QRM.Q());
+
+    // R1 (KP:p. 648)
+    // Find the first diagonal element index with (almost) zero value in R
+    label firstZeroElemi = 0;
+    {
+        const List<cmptType> diag(R.diag());
+
+        auto lessT = [&](const cmptType& x) { return mag(x) < tol; };
+
+        firstZeroElemi =
+            std::distance
+            (
+                diag.cbegin(),
+                std::find_if(diag.cbegin(), diag.cend(), lessT)
+            );
+    }
+
+    if (firstZeroElemi == 0)
+    {
+        FatalErrorInFunction
+            << "The largest (magnitude) diagonal element is (almost) zero."
+            << abort(FatalError);
+    }
+
+    // Exclude the first (almost) zero diagonal row and the rows below
+    // since R diagonal is already descending due to the QR column pivoting
+    const RectangularMatrix<cmptType> R1(R.subMatrix(0, 0, firstZeroElemi));
+
+    // R2 (KP:p. 648)
+    if (R1.n() < R1.m())
+    {
+        const SquareMatrix<cmptType> C(R1&R1);
+
+        QRMatrix<SquareMatrix<cmptType>> QRSolve
+        (
+            C,
+            QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+        );
+
+        RectangularMatrix<cmptType> R2
+        (
+            QRSolve.backSubstitution
+            (
+                QRSolve.R(),
+                QRSolve.Q() & (R1.T())
+            )
+        );
+
+        // R3 (KP:p. 648)
+        R2.resize(R.m(), R.n());
+
+        return MatrixType((QRM.P()^R2)^Q);
+    }
+    else
+    {
+        const SquareMatrix<cmptType> C(R1^R1);
+
+        QRMatrix<SquareMatrix<cmptType>> QRSolve
+        (
+            C,
+            QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+        );
+
+        RectangularMatrix<cmptType> R2
+        (
+            QRSolve.backSubstitution
+            (
+                QRSolve.R(),
+                QRSolve.Q() & R1
+            )
+        );
+
+        // R3
+        R2.resize(R.m(), R.n());
+
+        return MatrixType((QRM.P()^R2)^Q);
+    }
+}
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index 9058cd027e..415fbba853 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -28,19 +28,117 @@ Class
     Foam::QRMatrix
 
 Description
-    Class templated on matrix type to perform the QR decomposition using
-    Householder reflections on a square or rectangular matrix.
+    QRMatrix (i.e. QR decomposition, QR factorisation or orthogonal-triangular
+    decomposition) decomposes a scalar/complex matrix \c A into the following
+    matrix product:
+
+    \verbatim
+        A = Q*R,
+    \endverbatim
+
+    where
+     \c Q is a unitary similarity matrix,
+     \c R is an upper triangular matrix.
+
+    Reference:
+    \verbatim
+        QR decomposition:
+            mathworld.wolfram.com/QRDecomposition.html (Retrieved:17-06-19)
+            w.wiki/52X (Retrieved: 17-06-19)
+
+        QR decomposition with Householder reflector (tag:M):
+            Monahan, J. F. (2011).
+            Numerical methods of statistics.
+            Cambridge: Cambridge University Press.
+            DOI:10.1017/CBO9780511977176
+
+        QR decomposition with column pivoting (tag:QSB):
+            Quintana-Ortí, G., Sun, X., & Bischof, C. H. (1998).
+            A BLAS-3 version of the QR factorization with column pivoting.
+            SIAM Journal on Scientific Computing, 19(5), 1486-1494.
+            DOI:10.1137/S1064827595296732
+
+        Moore-Penrose inverse algorithm (tags:KP; KPP):
+            Katsikis, V. N., & Pappas, D. (2008).
+            Fast computing of the Moore-Penrose inverse matrix.
+            Electronic Journal of Linear Algebra, 17(1), 637-650.
+            DOI:10.13001/1081-3810.1287
+
+            Katsikis, V. N., Pappas, D., & Petralias, A. (2011).
+            An improved method for the computation of
+            the Moore–Penrose inverse matrix.
+            Applied Mathematics and Computation, 217(23), 9828-9834.
+            DOI:10.1016/j.amc.2011.04.080
+
+        Tolerance for the Moore-Penrose inverse algorithm (tag:TA):
+            Toutounian, F., & Ataei, A. (2009).
+            A new method for computing Moore–Penrose inverse matrices.
+            Journal of Computational and applied Mathematics, 228(1), 412-417.
+            DOI:10.1016/j.cam.2008.10.008
+
+        Operands of QR decomposition:
+            mathworld.wolfram.com/UnitaryMatrix.html (Retrieved:13-06-19)
+            mathworld.wolfram.com/UpperTriangularMatrix.html (Retrieved:16-06-19)
+            mathworld.wolfram.com/PermutationMatrix.html (Retrieved:20-06-19)
+
+        Back substitution:
+            mathworld.wolfram.com/GaussianElimination.html (Retrieved:15-06-19)
+    \endverbatim
+
+Usage
+    Input types:
+     - \c A can be a \c SquareMatrix<Type> or \c RectangularMatrix<Type>
+
+    Output types:
+     - \c Q is always of the type of the matrix \c A
+     - \c R is always of the type of the matrix \c A
+
+    Options for the output forms of \c QRMatrix (for an (m-by-n) input matrix
+    \c A with k = min(m, n)):
+     - outputTypes::FULL_R:     computes only \c R                   (m-by-n)
+     - outputTypes::FULL_QR:    computes both \c R and \c Q          (m-by-m)
+     - outputTypes::REDUCED_R:  computes only reduced \c R           (k-by-n)
+
+    Options where to store \c R:
+     - storeMethods::IN_PLACE:        replaces input matrix content with \c R
+     - storeMethods::OUT_OF_PLACE:    creates new object of \c R
+
+    Options for the computation of column pivoting:
+     - colPivoting::FALSE:            switches off column pivoting
+     - colPivoting::TRUE:             switches on column pivoting
+
+    Direct solution of linear systems A x = b is possible by solve() alongside
+    the following limitations:
+     - \c A         = a scalar square matrix
+     - output type  = outputTypes::FULL_QR
+     - store method = storeMethods::IN_PLACE
+
+Notes
+    - QR decomposition is not unique if \c R is not positive diagonal \c R.
+    - The option combination:
+      - outputTypes::REDUCED_R
+      - storeMethods::IN_PLACE
+      will not modify the rows of input matrix \c A after its nth row.
+    - Both FULL_R and REDUCED_R QR decompositions execute the same number of
+      operations. Yet REDUCED_R QR decomposition returns only the first n rows
+      of \c R if m > n for an input m-by-n matrix \c A.
+    - For m <= n, FULL_R and REDUCED_R will produce the same matrices.
+
+See also
+    Test-QRMatrix.C
 
 SourceFiles
-    QRMatrixI.H
     QRMatrix.C
+    QRMatrixI.H
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef QRMatrix_H
 #define QRMatrix_H
 
+#include "RectangularMatrix.H"
 #include "SquareMatrix.H"
+#include "complex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,32 +146,90 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class QRMatrix Declaration
+                            Class QRMatrix Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class MatrixType>
 class QRMatrix
 {
+
 public:
 
     typedef typename MatrixType::cmptType cmptType;
-    typedef SquareMatrix<cmptType> QMatrixType;
-    typedef MatrixType RMatrixType;
+    typedef SquareMatrix<cmptType> SMatrix;
+    typedef RectangularMatrix<cmptType> RMatrix;
+
+    //- Options for the output matrix forms of QRMatrix
+    enum outputTypes : uint8_t
+    {
+        FULL_R = 1,         //!< computes only \c R
+        FULL_QR = 2,        //!< computes both \c R and \c Q
+        REDUCED_R = 3       //!< computes only reduced \c R
+    };
+
+    //- Options where to store R
+    enum storeMethods : uint8_t
+    {
+        IN_PLACE = 1,       //!< replaces input matrix content with \c R
+        OUT_OF_PLACE = 2    //!< creates new object of \c R
+    };
+
+    //- Options for the computation of column pivoting
+    enum colPivoting : bool
+    {
+        FALSE = false,      //!< switches off column pivoting
+        TRUE  = true        //!< switches on column pivoting
+    };
 
 
 private:
 
     // Private Data
 
-        //- The Q-matrix
-        QMatrixType Q_;
+        //- Selected option for the output matrix forms of QRMatrix
+        outputTypes outputType_;
 
-        //- The R-matrix
-        RMatrixType R_;
+        //- Selected option where to store R
+        const storeMethods storeMethod_;
+
+        //- Selected option for the computation of column pivoting
+        const colPivoting colPivot_;
+
+        //- Unitary similarity matrix
+        MatrixType Q_;
+
+        //- Upper triangular matrix
+        MatrixType R_;
+
+        //- Permutation list (if column-pivoting is on)
+        labelList P_;
 
 
     // Private Member Functions
 
+        //- Compute the QR decomposition without the column pivoting
+        void qr
+        (
+            MatrixType& A
+        );
+
+        //- Compute the QR decomposition with the column pivoting
+        //  (QSB:Fig. 1)
+        void qrPivot
+        (
+            MatrixType& A
+        );
+
+        //- Compute output matrices in selected forms
+        //- using Householder reflectors
+        //  (M:Section 4)
+        inline void applyHouseholder
+        (
+            MatrixType& A,
+            const RMatrix& reflector,
+            const label k = 0
+        );
+
         //- Solve the linear system with the Field argument x initialized to
         //- the appropriate transformed source (e.g. Q.T()*source)
         //- and return the solution in x
@@ -90,27 +246,99 @@ private:
         ) const;
 
 
+protected:
+
+    // Protected Member Functions
+
+        //- Compute Householder reflector on a given matrix column, u
+        //  (M:Eq.6.814)
+        inline RMatrix householderReflector
+        (
+            RMatrix u
+        );
+
+        //- Apply (in-place) Householder reflectors from the left side: u*A
+        inline void applyLeftReflector
+        (
+            MatrixType& A,
+            const RMatrix& reflector,
+            const label k = 0,
+            const label k1 = 0
+        );
+
+        //- Apply (in-place) Householder reflectors from the right side: (u*A)*u
+        inline void applyRightReflector
+        (
+            MatrixType& A,
+            const RMatrix& reflector,
+            const label k = 0
+        );
+
+
 public:
 
     // Constructors
 
         //- Construct null
-        inline QRMatrix();
+        QRMatrix();
 
-        //- Construct decomposing given matrix
-        QRMatrix(const MatrixType& M);
+        //- Construct QRMatrix without performing the decomposition
+        QRMatrix
+        (
+            const outputTypes outputType,
+            const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
+            const colPivoting colPivot = colPivoting::FALSE
+        );
+
+        //- Construct QRMatrix and perform the QR decomposition
+        QRMatrix
+        (
+            MatrixType& A,
+            const outputTypes outputType,
+            const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
+            const colPivoting colPivot = colPivoting::FALSE
+        );
+
+        //- Construct QRMatrix and perform the QR decomposition
+        QRMatrix
+        (
+            const MatrixType& A,
+            const outputTypes outputType,
+            const colPivoting colPivot = colPivoting::FALSE
+        );
 
 
     // Member Functions
 
-        //- Return Q-matrix
-        inline const QMatrixType& Q() const;
+        // Information
 
-        //- Return R-matrix
-        inline const RMatrixType& R() const;
+        //- Return the unitary similarity matrix
+        //  Includes implicit round-to-zero as mutable operation
+        inline const MatrixType& Q() const;
 
-        //- Decompose given matrix
-        void decompose(const MatrixType& M);
+        //- Return the upper triangular matrix
+        inline const MatrixType& R() const;
+
+        //- Return the permutation order (P) as a list
+        inline const labelList& orderP() const;
+
+        //- Create and return the permutation matrix
+        inline SMatrix P() const;
+
+
+        // Algorithm
+
+        //- Compute QR decomposition according to constructor settings
+        void decompose
+        (
+            MatrixType& A
+        );
+
+        //- Compute QR decomposition according to constructor settings
+        void decompose
+        (
+            const MatrixType& A
+        );
 
         //- Solve the linear system with the given source
         //- and return the solution in the argument x
@@ -145,10 +373,38 @@ public:
         ) const;
 
         //- Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source
-        QMatrixType inv() const;
+        SMatrix inv() const;
+
+        //- Solve a row-echelon-form linear system starting from the bottom
+        //  Back substitution: Ax = b where:
+        //  A = Non-singular upper-triangular square matrix (m-by-m)
+        //  x = Solution (m-by-p)
+        //  b = Source (m-by-p)
+        RMatrix backSubstitution
+        (
+            const SMatrix& A,
+            const RMatrix& b
+        );
 };
 
 
+// * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * * //
+
+//- (Out-of-place) Moore-Penrose inverse of singular/non-singular
+//- square/rectangular scalar/complex matrices
+//  (KPP:p. 9834; KP:p. 648)
+//  The tolerance (i.e. tol) to ensure the R1 matrix full-rank is given as 1e-13
+//  in the original paper (KPP:p. 9834).  Nonetheless, the tolerance = 1e-5
+//  given by (TA; mentioned in (KPP:p. 9832)) was observed to be more robust
+//  for the tested scenarios.
+template<class MatrixType>
+MatrixType pinv
+(
+    const MatrixType& A,
+    const scalar tol = 1e-5
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
index 69102fa722..929f8b6fc8 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
@@ -26,29 +26,167 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class MatrixType>
-inline Foam::QRMatrix<MatrixType>::QRMatrix()
-{}
+inline void Foam::QRMatrix<MatrixType>::applyHouseholder
+(
+    MatrixType& A,
+    const RMatrix& reflector,
+    const label k
+)
+{
+    applyLeftReflector(A, reflector, k, k);
+
+    if (outputType_ == outputTypes::FULL_QR)
+    {
+        applyRightReflector(Q_, reflector, k);
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+template<class MatrixType>
+inline Foam::RectangularMatrix<typename MatrixType::cmptType>
+Foam::QRMatrix<MatrixType>::householderReflector
+(
+    RMatrix u
+)
+{
+    #ifdef FULLDEBUG
+    // Check if the given RectangularMatrix is effectively a column vector
+    if (u.n() != 1)
+    {
+        FatalErrorInFunction
+            << "Input matrix is not a column vector." << exit(FatalError);
+    }
+    #endif
+
+    scalar magnitude(mag(u(0,0)));
+    if (magnitude < VSMALL)
+    {
+        magnitude = SMALL;
+        #if FULLDEBUG
+        FatalErrorInFunction
+            << "Almost zero leading elem in Householder matrix."
+            << abort(FatalError);
+        #endif
+    }
+
+    u(0,0) += u(0,0)/magnitude*u.columnNorm(0);
+
+    scalar colNorm(u.columnNorm(0));
+    if (colNorm < VSMALL)
+    {
+        colNorm = SMALL;
+        #if FULLDEBUG
+        FatalErrorInFunction
+            << "Almost zero norm in the Householder matrix."
+            << abort(FatalError);
+        #endif
+    }
+
+    u /= cmptType(colNorm);
+
+    return u;
+}
+
+
+template<class MatrixType>
+inline void Foam::QRMatrix<MatrixType>::applyLeftReflector
+(
+    MatrixType& A,
+    const RMatrix& reflector,
+    const label k,
+    const label k1
+)
+{
+    //  const RMatrix& A0(A.subMatrix(k1, k));
+    //  A0 -= (cmptType(2)*reflector)*(reflector & A0);
+
+    for (label j = k; j < A.n(); ++j)
+    {
+        cmptType sum = Zero;
+        for (label i = 0; i < reflector.m(); ++i)
+        {
+            sum += Detail::conj(reflector(i, 0))*A(i + k1, j);
+        }
+
+        sum *= cmptType(2);
+        for (label i = 0; i < reflector.m(); ++i)
+        {
+            A(i + k1, j) -= reflector(i, 0)*sum;
+        }
+    }
+}
+
+
+template<class MatrixType>
+inline void Foam::QRMatrix<MatrixType>::applyRightReflector
+(
+    MatrixType& A,
+    const RMatrix& reflector,
+    const label k
+)
+{
+    // const RMatrix A0(A.subMatrix(0, k));
+    // A0 -= ((A0*reflector)^(cmptType(2)*reflector))
+
+    for (label i = 0; i < A.m(); ++i)
+    {
+        cmptType sum = Zero;
+        for (label j = 0; j < reflector.m(); ++j)
+        {
+            sum += A(i, j + k)*reflector(j, 0);
+        }
+
+        sum *= cmptType(2);
+        for (label j = 0; j < reflector.m(); ++j)
+        {
+            A(i, j + k) -= Detail::conj(reflector(j, 0))*sum;
+        }
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MatrixType>
-inline const typename Foam::QRMatrix<MatrixType>::QMatrixType&
-Foam::QRMatrix<MatrixType>::Q() const
+inline const MatrixType& Foam::QRMatrix<MatrixType>::Q() const
 {
+    const_cast<MatrixType&>(Q_).round();
+
     return Q_;
 }
 
 
 template<class MatrixType>
-inline const typename Foam::QRMatrix<MatrixType>::RMatrixType&
-Foam::QRMatrix<MatrixType>::R() const
+inline const MatrixType& Foam::QRMatrix<MatrixType>::R() const
 {
     return R_;
 }
 
 
+template<class MatrixType>
+inline const Foam::labelList& Foam::QRMatrix<MatrixType>::orderP() const
+{
+    return P_;
+}
+
+
+template<class MatrixType>
+inline Foam::SquareMatrix<typename Foam::QRMatrix<MatrixType>::cmptType>
+Foam::QRMatrix<MatrixType>::P() const
+{
+    SquareMatrix<cmptType> permMat(P_.size(), Zero);
+
+    forAll(P_, jcol)
+    {
+        permMat(P_[jcol], jcol) = pTraits<cmptType>::one;
+    }
+
+    return permMat;
+}
+
+
 // ************************************************************************* //