From 1ef855cf0bb2dac72aa4701276ee68ddacf09e82 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Fri, 20 May 2022 15:08:43 +0100
Subject: [PATCH] ENH: QRMatrix: refactor the QR decomposition algorithms
 (fixes #2212)

---
 applications/test/TestTools/TestTools.H       |   91 +-
 .../test/matrices/QRMatrix/Make/options       |    2 +-
 .../test/matrices/QRMatrix/Test-QRMatrix.C    | 1613 ++++-------------
 src/OpenFOAM/matrices/QRMatrix/QRMatrix.C     |  633 ++++---
 src/OpenFOAM/matrices/QRMatrix/QRMatrix.H     |  350 ++--
 src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H    |  158 +-
 6 files changed, 980 insertions(+), 1867 deletions(-)

diff --git a/applications/test/TestTools/TestTools.H b/applications/test/TestTools/TestTools.H
index 3fed3b7d5d..82dfe5e38b 100644
--- a/applications/test/TestTools/TestTools.H
+++ b/applications/test/TestTools/TestTools.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,7 @@ using namespace Foam;
 #include "complex.H"
 #include "Matrix.H"
 #include "Random.H"
+#include <chrono>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,6 +47,27 @@ unsigned nTest_ = 0;
 unsigned nFail_ = 0;
 
 
+// Return the absolute tolerance value for bitwise comparisons of floatScalars
+floatScalar getTol(floatScalar)
+{
+    return 1e-2;
+}
+
+
+// Return the absolute tolerance value for bitwise comparisons of doubleScalars
+doubleScalar getTol(doubleScalar)
+{
+    return 1e-10;
+}
+
+
+// Return the absolute tolerance value for bitwise comparisons of doubleScalars
+doubleScalar getTol(complex)
+{
+    return 1e-10;
+}
+
+
 // Create a non-complex random Matrix.
 template<class MatrixType>
 typename std::enable_if
@@ -174,10 +196,14 @@ typename std::enable_if
     const Type& x,
     const Type& y,
     const scalar absTol = 0,    //<! useful for cmps near zero
-    const scalar relTol = 1e-8  //<! are values the same within 8 decimals
+    const scalar relTol = 1e-8, //<! are values the same within 8 decimals
+    const bool verbose = false
 )
 {
-    Info<< msg << x << "?=" << y << endl;
+    if (verbose)
+    {
+        Info<< msg << x << "?=" << y << endl;
+    }
 
     unsigned nFail = 0;
 
@@ -212,10 +238,14 @@ typename std::enable_if
     const Type& x,
     const Type& y,
     const scalar absTol = 0,
-    const scalar relTol = 1e-8
+    const scalar relTol = 1e-8,
+    const bool verbose = false
 )
 {
-    Info<< msg << x << "?=" << y << endl;
+    if (verbose)
+    {
+        Info<< msg << x << "?=" << y << endl;
+    }
 
     unsigned nFail = 0;
 
@@ -253,10 +283,14 @@ typename std::enable_if
     const Type1& x,
     const Type2& y,
     const scalar absTol = 0,
-    const scalar relTol = 1e-8
+    const scalar relTol = 1e-8,
+    const bool verbose = false
 )
 {
-    Info<< msg << x << "?=" << y << endl;
+    if (verbose)
+    {
+        Info<< msg << x << "?=" << y << endl;
+    }
 
     unsigned nFail = 0;
 
@@ -284,10 +318,14 @@ void cmp
 (
     const word& msg,
     const bool x,
-    const bool y
+    const bool y,
+    const bool verbose = false
 )
 {
-    Info<< msg << x << "?=" << y << endl;
+    if (verbose)
+    {
+        Info<< msg << x << "?=" << y << endl;
+    }
 
     unsigned nFail = 0;
 
@@ -306,4 +344,39 @@ void cmp
 }
 
 
+// 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;
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/test/matrices/QRMatrix/Make/options b/applications/test/matrices/QRMatrix/Make/options
index 4e772fdf9d..b23590707c 100644
--- a/applications/test/matrices/QRMatrix/Make/options
+++ b/applications/test/matrices/QRMatrix/Make/options
@@ -1,2 +1,2 @@
-/* EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude */
+EXE_INC = -I../../TestTools
 /* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/QRMatrix/Test-QRMatrix.C b/applications/test/matrices/QRMatrix/Test-QRMatrix.C
index 6790727152..9e51d67eb0 100644
--- a/applications/test/matrices/QRMatrix/Test-QRMatrix.C
+++ b/applications/test/matrices/QRMatrix/Test-QRMatrix.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -23,946 +23,338 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
+Application
+    Test-QRMatrix
+
+Description
+    Tests for \c QRMatrix constructors, and member functions
+    using \c doubleScalar and \c complex base types.
+
 \*---------------------------------------------------------------------------*/
 
 #include "MatrixTools.H"
 #include "QRMatrix.H"
-#include "Random.H"
+#include "complex.H"
 #include "IOmanip.H"
+#include "TestTools.H"
+
+#define PRINTALL false
 
 using namespace Foam;
-using namespace Foam::MatrixTools;
-
-#define equal MatrixTools::equal
-#define RUNALL true
-const bool verbose = true;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void horizontalLine()
+// Verify if Q is a unitary matrix for rectangular matrices
+template<class Type>
+void verify_Q(const RectangularMatrix<Type>& Q)
 {
-    Info<< "+---------+---------+---------+---------+---------+" << nl;
+    // non-square unitary matrices are not possible - do nothing
 }
 
 
-// Create random scalar-type matrix
-template<class MatrixType>
-typename std::enable_if
-<
-    !std::is_same<complex, typename MatrixType::cmptType>::value,
-    MatrixType
->::type makeRandomMatrix
+// Verify if Q is a unitary matrix for square matrices
+template<class Type>
+void verify_Q
 (
-    const labelPair& dims,
-    Random& rndGen
+    const SquareMatrix<Type>& Q
 )
 {
-    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);
+    const SquareMatrix<Type> QTQ(Q&Q);
+    const SquareMatrix<Type> QQT(Q^Q);
+    SquareMatrix<Type> IMatrix(Q.sizes(), Zero);
     for (label i = 0; i < IMatrix.n(); ++i)
     {
-        IMatrix(i,i) = pTraits<typename MatrixType::cmptType>::one;
+        IMatrix(i,i) = pTraits<Type>::one;
     }
-    equal(QTQ, IMatrix, verbose);
-    equal(QQT, IMatrix, verbose);
-    equal(QTQ, QQT, verbose);
+    cmp("# Q.T()*Q = I: ", flt(QTQ), flt(IMatrix), getTol(Type(0)));
+    cmp("# Q*Q.T() = I: ", flt(QQT), flt(IMatrix), getTol(Type(0)));
+    cmp("# QTQ = QQT: ", flt(QTQ), flt(QQT), getTol(Type(0)));
 }
 
 
-// Checks if R matrix is an upper-triangular matrix
+// Verify if R is an upper-triangular matrix
 template<class MatrixType>
-void cross_check_QRMatrix
+void verify_R
 (
     const MatrixType& R
 )
 {
-    InfoNumPyFormat(R, "R");
-
+    DynamicList<scalar> ltriR;
+    label k = 0;
     // 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);
+                ltriR.append(mag(R(i, j)));
+                ++k;
             }
         }
     }
+    const List<scalar> ltri(k, Zero);
+    cmp("# R(i, j) = 0 for i > j: ", ltri, ltriR);
 }
 
 
-// Checks if the given matrix can be reconstructed by Q and R matrices
+// Verify if R is an upper-triangular matrix with diag elems non-increasing
 template<class MatrixType>
-void cross_check_QRMatrix
+void verify_R_pivoting
 (
-    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);
-    }
+    cmp
+    (
+        "# Diag elems of R non-increasing: |R11| >= |R22| >= ... >= |Rnn|: ",
+        diag0,
+        diag1
+    );
 }
 
 
-// Checks if the given matrix can be reconstructed by column-pivoted Q and R
+// Verify if the given matrix can be reconstructed by Q and R matrices
 template<class MatrixType>
-void cross_check_QRMatrix
+void verify_QR
 (
-    const MatrixType& Aorig,
+    const MatrixType& A,
     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;
+    // mathworld.wolfram.com/QRDecomposition.html (Retrieved:18-06-19)
     const MatrixType AReconstruct(Q*R);
-    const MatrixType AP(Aorig*P);
-    equal(AP, AReconstruct, verbose);
+    cmp("# Q*R = A: ", flt(A), flt(AReconstruct));
 }
 
 
-// Checks each constructor of QRMatrix type
+// Verify if the given matrix can be reconstructed by Q, R and P matrices
 template<class MatrixType>
-void verification_QRMatrix
+void verify_QR_pivoting
 (
-    MatrixType& A
+    const MatrixType& A,
+    const MatrixType& Q,
+    const MatrixType& R,
+    const SquareMatrix<typename MatrixType::cmptType>& P
 )
 {
-    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
+    // w.wiki/54m (Retrieved:20-06-19)
+    const MatrixType AReconstruct(Q*R);
+    const MatrixType AP(A*P);
+    cmp("# Q*R = A*P: ", flt(AReconstruct), flt(AP));
 }
 
 
-// Checks the direct solution of a linear system by QRMatrix
+// Verify the given QR decomposition
 template<class MatrixType>
-void verification_QRMatrixSolve
+void verify
 (
-    MatrixType& A
+    const QRMatrix<MatrixType>& QRM,
+    const MatrixType& A
 )
 {
-    typedef SquareMatrix<typename MatrixType::cmptType> SMatrix;
+    const MatrixType& Q = QRM.Q();
+    const MatrixType& R = QRM.R();
 
-    // Create working copies of matrix A
-    const MatrixType Aorig(A);
+    #if PRINTALL
+    InfoNumPyFormat(A, "A");
+    InfoNumPyFormat(Q, "Q");
+    InfoNumPyFormat(R, "R");
+    #endif
+
+    verify_Q(Q);
+    verify_R(R);
+    verify_QR(A, Q, R);
+    Info<< endl;
+}
+
+
+// Verify the given QR decomposition with column pivoting
+template<class MatrixType>
+void verify_pivoting
+(
+    const QRMatrix<MatrixType>& QRM,
+    const MatrixType& A
+)
+{
+    const MatrixType& Q = QRM.Q();
+    const MatrixType& R = QRM.R();
+    const SquareMatrix<typename MatrixType::cmptType>& P = QRM.P();
+
+    #if PRINTALL
+    InfoNumPyFormat(A, "A");
+    InfoNumPyFormat(Q, "Q");
+    InfoNumPyFormat(R, "R");
+    InfoNumPyFormat(P, "P");
+    #endif
+
+    verify_Q(Q);
+    verify_R(R);
+    verify_R_pivoting(R);
+    verify_QR_pivoting(A, Q, R, P);
+    Info<< endl;
+}
+
+
+// Verify various QR decompositions
+template<class MatrixType>
+void verify_decomposition
+(
+    const MatrixType& A
+)
+{
+    Info<< "## Verify various QR decompositions" << nl << endl;
+    {
+        Info<< "# modes::FULL" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            A,
+            QRMatrix<MatrixType>::modes::FULL
+        );
+        verify(QRM, A);
+    }
+    {
+        Info<< "# modes::ECONOMY" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            A,
+            QRMatrix<MatrixType>::modes::ECONOMY
+        );
+        verify(QRM, A);
+    }
+    {
+        Info<< "# modes::FULL, outputs::BOTH_QR, column-pivoting" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            A,
+            QRMatrix<MatrixType>::modes::FULL,
+            QRMatrix<MatrixType>::outputs::BOTH_QR,
+            QRMatrix<MatrixType>::pivoting::TRUE
+        );
+        verify_pivoting(QRM, A);
+    }
+    {
+        Info<< "# modes::ECONOMY, outputs::BOTH_QR, column-pivoting" << nl;
+        QRMatrix<MatrixType> QRM
+        (
+            A,
+            QRMatrix<MatrixType>::modes::ECONOMY,
+            QRMatrix<MatrixType>::outputs::BOTH_QR,
+            QRMatrix<MatrixType>::pivoting::TRUE
+        );
+        verify_pivoting(QRM, A);
+    }
+}
+
+
+// Verify the Moore-Penrose pseudo-inverse function
+template<class MatrixType>
+void verify_pinv
+(
+    const MatrixType& A
+)
+{
+    const scalar absTol = 1e-6;
+
+    Info<< "## Verify Moore-Penrose pseudo-inverse: pinv()" << nl << endl;
+    const MatrixType APlus(MatrixTools::pinv(A));
+
+    #if PRINTALL
+    InfoNumPyFormat(A, "A");
+    InfoNumPyFormat(APlus, "A^+");
+    #endif
+
+    {
+        const MatrixType APlusPlus(MatrixTools::pinv(APlus));
+        cmp("# (A^+)^+ = A: ", flt(APlusPlus), flt(A), absTol);
+    }
+
+    // The Moore-Penrose conditions
+    // w.wiki/5RL (Retrieved: 29-06-19)
+    {
+        const MatrixType AAPlusA((A*APlus)*A);
+        cmp("# A*(A^+)*A = A: ", flt(AAPlusA), flt(A), absTol);
+    }
+
+    {
+        const MatrixType APlusAAPlus((APlus*A)*APlus);
+        cmp("# (A^+)*A*(A^+) = A: ", flt(APlusAAPlus), flt(APlus), absTol);
+    }
+
+    {
+        const MatrixType AAPlus(A*APlus);
+        const MatrixType AAPlusT(AAPlus.T());
+        cmp("# (A*(A^+)).T() = A*(A^+): ", flt(AAPlusT), flt(AAPlus), absTol);
+    }
+
+    {
+        const MatrixType APlusA(APlus*A);
+        const MatrixType APlusAT(APlusA.T());
+        cmp("# ((A^+)*A = ((A^+)*A).T(): ", flt(APlusA), flt(APlusAT), absTol);
+    }
+}
+
+
+// Verify the direct solution of a linear system by QRMatrix
+void verify_solve
+(
+    const SquareMatrix<complex>& A
+)
+{
+    // do nothing
+}
+
+
+// Verify the direct solution of a linear system by QRMatrix
+void verify_solve
+(
+    const SquareMatrix<scalar>& A
+)
+{
+    Info<< "## Verify direct solution of A*x = b with A = Q*R:" << nl << endl;
+
+    typedef SquareMatrix<scalar> SMatrix;
+    const scalar absTol = 1e-6;
 
     // Direct solution of a linear system by QR decomposition
-    #if (0 | RUNALL)
-    {
-        MatrixType A0(A);
-        QRMatrix<MatrixType> QRM
-        (
-            A0,
-            QRMatrix<MatrixType>::outputTypes::FULL_QR
-        );
+    QRMatrix<SMatrix> QRM
+    (
+        A,
+        QRMatrix<SMatrix>::modes::FULL,
+        QRMatrix<SMatrix>::outputs::BOTH_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 residual(A.m(), 0);
+    const scalarField source(A.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);
+    const scalarField x(QRM.solve(source));
+    const scalarField AxMinusSource(A*x - source);
+    const scalarField xMinusQRinvSource(x - QRM.inv()*source);
+    const SMatrix QRinvA(QRM.inv()*A);
+    const SMatrix identity(A.m(), I);
 
-        forAll(residual, i)
-        {
-            isEqual(A0xMinusSource[i], residual[i], verbose);
-            isEqual(xMinusQRinvSource[i], residual[i], verbose);
-        }
-        equal(QRinvA0, identity, verbose);
-    }
-    #endif
+    cmp("# A*x - b = residual: ", AxMinusSource, residual, absTol);
+    cmp("# x - inv(QR)*b = residual: ", xMinusQRinvSource, residual, absTol);
+    cmp("# inv(QR)*A = I: ", flt(QRinvA), flt(identity), absTol);
 }
 
 
-// Checks the parallel tall-skinny QR decomposition
+// Verify the parallel tall-skinny QR decomposition
 template<class Type>
-void verification_tsqr
+void verify_tsqr
 (
     Random& rndGen
 )
 {
+    Info<< "## Verify parallel direct tall-skinny QR decomposition:"
+        << nl << endl;
+
     typedef RectangularMatrix<Type> RMatrix;
 
     // Size of the full matrix and its partitions
@@ -978,7 +370,7 @@ void verification_tsqr
         mRowsSum += mRows;
     }
 
-    # if 0
+    #if PRINTALL
     Info<< "mRowsSum = " << mRowsSum << nl;
     Info<< "nColsSum = " << nColsSum << nl;
     Info<< "mRowsList = " << mRowsList << nl;
@@ -999,8 +391,9 @@ void verification_tsqr
     QRMatrix<RMatrix> QRM
     (
         A,
-        QRMatrix<RMatrix>::outputTypes::REDUCED_R,
-        QRMatrix<RMatrix>::colPivoting::FALSE
+        QRMatrix<RMatrix>::modes::ECONOMY,
+        QRMatrix<RMatrix>::outputs::ONLY_R,
+        QRMatrix<RMatrix>::pivoting::FALSE
     );
     const RMatrix& R = QRM.R();
 
@@ -1023,8 +416,9 @@ void verification_tsqr
         QRMatrix<RMatrix> subQRM
         (
             subA,
-            QRMatrix<RMatrix>::outputTypes::REDUCED_R,
-            QRMatrix<RMatrix>::colPivoting::FALSE
+            QRMatrix<RMatrix>::modes::ECONOMY,
+            QRMatrix<RMatrix>::outputs::ONLY_R,
+            QRMatrix<RMatrix>::pivoting::FALSE
         );
         const RMatrix& subR = subQRM.R();
 
@@ -1036,8 +430,9 @@ void verification_tsqr
     QRMatrix<RMatrix> TSQR
     (
         subRs,
-        QRMatrix<RMatrix>::outputTypes::REDUCED_R,
-        QRMatrix<RMatrix>::colPivoting::FALSE
+        QRMatrix<RMatrix>::modes::ECONOMY,
+        QRMatrix<RMatrix>::outputs::ONLY_R,
+        QRMatrix<RMatrix>::pivoting::FALSE
     );
     const RMatrix& TSQRR = TSQR.R();
 
@@ -1048,433 +443,161 @@ void verification_tsqr
     }
 
     // Compare magnitude of R, since R is not unique up to the sign
-    equal(RMag, TSQRRMag, verbose);
+    cmp("# TSQR: ", flt(RMag), flt(TSQRRMag));
 
-    #if 0
+    #if PRINTALL
     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
-)
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+void test_constructors(MatrixType)
 {
-    // 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");
+        Info<< "# Construct from a Matrix and modes" << nl;
+        MatrixType A({5,5}, Zero);
+        const QRMatrix<MatrixType> QRM
+        (
+            A,
+            QRMatrix<MatrixType>::modes::FULL
+        );
     }
-    #endif
-
-    #if (0 | RUNALL)
     {
-        Info<< nl << "# A*X = b:" << nl;
-        const RMatrix AX(A*X);
-        equal(AX, b, verbose, 10, 1e-3, 1e-3);
+        Info<< "# Construct from a const Matrix and modes" << nl;
+        const MatrixType A({5,5}, Zero);
+        const QRMatrix<MatrixType> QRM
+        (
+            A,
+            QRMatrix<MatrixType>::modes::FULL
+        );
     }
-    #endif
 }
 
 
-// Checks the Moore-Penrose pseudo-inverse function
 template<class MatrixType>
-void verification_pinv
-(
-    const MatrixType& A
-)
+void test_decomposition(MatrixType)
 {
-    Info<< "### Moore-Penrose pseudo-inverse: pinv()" << nl << nl;
-    const MatrixType APlus(pinv(A));
+    typedef typename MatrixType::cmptType cmptType;
+    typedef SquareMatrix<cmptType> SMatrix;
+    typedef RectangularMatrix<cmptType> RMatrix;
 
-    #if 0
-    InfoNumPyFormat(A, "A");
-    InfoNumPyFormat(APlus, "A^+");
-    #endif
+    Random rndGen(1234);
+    const label szTests = 20;
+    const label min = 10;
+    const label max = 50;
 
-    if (0 | RUNALL)
     {
-        Info<< nl << "# (A^+)^+ = A:" << nl;
-        const MatrixType APlusPlus(pinv(APlus));
-        equal(APlusPlus, A, verbose);
+        for (label i = 0; i < szTests; ++i)
+        {
+            const label m = rndGen.position(min, max);
+            const label n = rndGen.position(min, max);
+
+            const RMatrix A
+            (
+                makeRandomMatrix<RMatrix>({m,n}, rndGen)
+            );
+
+            verify_decomposition(A);
+            verify_pinv(A);
+        }
     }
 
-    // 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);
+        for (label i = 0; i < szTests; ++i)
+        {
+            const label m = rndGen.position(min, max);
+
+            const SMatrix A
+            (
+                makeRandomMatrix<SMatrix>({m,m}, rndGen)
+            );
+
+            verify_decomposition(A);
+            verify_pinv(A);
+            verify_solve(A);
+        }
     }
 
-    if (0 | RUNALL)
+    for (label i = 0; i < szTests; ++i)
     {
-        Info<< nl << "# (A^+)*A*(A^+) = A:" << nl;
-        const MatrixType APlusAAPlus((APlus*A)*APlus);
-        equal(APlusAAPlus, APlus, verbose);
+        verify_tsqr<cmptType>(rndGen);
     }
+}
 
-    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);
-    }
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test decomposition: "<< typeID[I] <<" ##" << nl;
+    test_decomposition(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
 }
 
 
 // * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
 
-int main(int argc, char *argv[])
+int main()
 {
-    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;
+    const std::tuple
+    <
+        RectangularMatrix<doubleScalar>,
+        RectangularMatrix<complex>,
+        SquareMatrix<doubleScalar>,
+        SquareMatrix<complex>
+    > types
+    (
+        std::make_tuple(Zero, Zero, Zero, Zero)
+    );
 
-    // SquareMatrix<scalar>
-    #if (0 | RUNALL)
+    const List<word> typeID
+    ({
+        "RectangularMatrix<doubleScalar>",
+        "RectangularMatrix<complex>",
+        "SquareMatrix<doubleScalar>",
+        "SquareMatrix<complex>"
+    });
+
+    std::chrono::steady_clock::time_point begin =
+        std::chrono::steady_clock::now();
+
+    run_tests(types, typeID);
+
+    std::chrono::steady_clock::time_point end =
+        std::chrono::steady_clock::now();
+    Info<< "Time elapsed = "
+        << std::chrono::duration_cast<std::chrono::microseconds>(end - begin)
+            .count()
+        << "[µs]" << endl;
+
+    if (nFail_)
     {
-        horizontalLine();
-
-        Info<< "## " << numberOfTests << " SquareMatrix<scalar> A:" << nl;
-
-        for (label i = 0; i < numberOfTests; ++i)
-        {
-            const label mRows = rndGen.position<label>(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();
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
     }
-    #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<label>(1, 50);
-            const label nCols = rndGen.position<label>(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<label>(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<label>(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<label>(1, 50);
-            const label nCols = rndGen.position<label>(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<label>(20, 50);
-                const label pCols = rndGen.position<label>(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<label>(20, 50);
-                const label pCols = rndGen.position<label>(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;
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
 
     return 0;
 }
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index 0b63f2017b..f7bf7db3eb 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,96 +31,292 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::qr
+Foam::label Foam::QRMatrix<MatrixType>::calcShapeFactor
 (
-    MatrixType& A
-)
+    const MatrixType& A
+) const
 {
-    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)
+    if (mode_ == modes::FULL)
     {
-        const RMatrix u(householderReflector(A.subColumn(k, k)));
-
-        applyHouseholder(A, u, k);
+        return A.m();
+    }
+    else if (mode_ == modes::ECONOMY)
+    {
+        return min(A.m(), A.n());
     }
 
-    if (outputType_ == outputTypes::REDUCED_R)
-    {
-        A.resize(A.n(), A.n());
-    }
+    return 0;
 }
 
 
 template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::qrPivot
+void Foam::QRMatrix<MatrixType>::decompose
 (
-    MatrixType& A
+    MatrixType& AT
 )
 {
-    const label nCols = A.n();
-    const label nIter = min(A.m() - 1, A.n());
+    const label n = AT.m();
+    const label m = AT.n();
 
-    // Initialise permutation vector, and column norm vector
-    P_ = identity(nCols);
+    List<cmptType> Rdiag(n, Zero);
 
-    // Initialise vector norms of each column of A
-    List<scalar> colNorms(nCols);
-    for (label k = 0; k < nCols; ++k)
+    for (label k = 0; k < n; ++k)
     {
-        colNorms[k] = A.columnNorm(k, true);
+        // Compute 2-norm of k-th column without under/overflow
+        scalar nrm = 0;
+
+        for (label i = k; i < m; ++i)
+        {
+            nrm = Foam::hypot(nrm, mag(AT(k,i)));
+        }
+
+        // Avoid divide-by-zero. Use compare with == 0 not mag or VSMALL etc,
+        // otherwise wrong results (may need more investigation)
+        if (nrm != scalar(0))
+        {
+            // Form k-th Householder vector
+            if (mag(AT(k,k)) < 0)
+            {
+                nrm *= -1;
+            }
+
+            for (label i = k; i < m; ++i)
+            {
+                AT(k,i) /= nrm;
+            }
+
+            AT(k,k) += scalar(1);
+
+            // Apply transformation to remaining columns
+            for (label j = k + 1; j < n; ++j)
+            {
+                cmptType s(Zero);
+
+                for (label i = k; i < m; ++i)
+                {
+                    s += Detail::conj(AT(k,i))*AT(j,i);
+                }
+
+                if (mag(AT(k,k)) > SMALL)
+                {
+                    s /= -AT(k,k);
+                }
+
+                for (label i = k; i < m; ++i)
+                {
+                    AT(j,i) += s*AT(k,i);
+                }
+            }
+        }
+
+        Rdiag[k] = -nrm;
     }
 
-    // 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);
+    calcQ(AT);
 
-        // Column pivoting
-        const label maxColNormi =
+    calcR(AT, Rdiag);
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose
+(
+    MatrixType& AT,
+    const bool pivoting
+)
+{
+    const label n = AT.m();
+    const label m = AT.n();
+    const label sz = min(m,n);
+
+    // Initialise permutation vector, and column norm vector
+    p_ = identity(n);
+
+    // Initialise vector norms of each column of A
+    List<scalar> norms(n, Zero);
+    for (label k = 0; k < n; ++k)
+    {
+        for (label i = 0; i < m; ++i)
+        {
+            norms[k] += magSqr(AT(k, i));
+        }
+    }
+
+    List<cmptType> Rdiag(n, Zero);
+
+    for (label k = 0; k < sz; ++k)
+    {
+        const auto it = std::next(norms.cbegin(), k);
+        const label maxNormi =
             std::distance
             (
-                subColNorms.cbegin(),
-                std::max_element(subColNorms.cbegin(), subColNorms.cend())
+                it,
+                std::max_element(it, norms.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)))
+        // Swap A, P and norms according to pivot column if the current
+        // column is not the max norm column. Also avoid any column swaps
+        // where the leading element is very small
+        if (mag(AT(k + maxNormi,k)) > SMALL && maxNormi != 0)
         {
-            const RMatrix R1(A.subColumn(k));
-            const RMatrix R2(A.subColumn(maxColNormi + k));
-            A.subColumn(k) = R2;
-            A.subColumn(maxColNormi + k) = R1;
+            const RMatrix R1(AT.subRow(k));
+            const RMatrix R2(AT.subRow(maxNormi + k));
 
-            Swap(P_[k], P_[maxColNormi + k]);
-            Swap(colNorms[k], colNorms[maxColNormi + k]);
+            AT.subRow(k) = R2;
+            AT.subRow(maxNormi + k) = R1;
+
+            Swap(p_[k], p_[maxNormi + k]);
+            Swap(norms[k], norms[maxNormi + k]);
         }
 
         {
-            const RMatrix u(householderReflector(A.subColumn(k, k)));
+            // Compute 2-norm of k-th column without under/overflow
+            scalar nrm = 0;
 
-            applyHouseholder(A, u, k);
+            for (label i = k; i < m; ++i)
+            {
+                nrm = Foam::hypot(nrm, mag(AT(k,i)));
+            }
+
+            if (nrm != scalar(0))
+            {
+                // Form k-th Householder vector
+                if (mag(AT(k,k)) < 0)
+                {
+                    nrm *= -1;
+                }
+
+                for (label i = k; i < m; ++i)
+                {
+                    AT(k,i) /= nrm;
+                }
+
+                AT(k,k) += scalar(1);
+
+                // Apply transformation to remaining columns
+                for (label j = k + 1; j < n; ++j)
+                {
+                    cmptType s(Zero);
+
+                    for (label i = k; i < m; ++i)
+                    {
+                        s += Detail::conj(AT(k,i))*AT(j,i);
+                    }
+
+                    if (mag(AT(k,k)) > SMALL)
+                    {
+                        s /= -AT(k,k);
+                    }
+
+                    for (label i = k; i < m; ++i)
+                    {
+                        AT(j,i) += s*AT(k,i);
+                    }
+                }
+            }
+
+            Rdiag[k] = -nrm;
         }
 
         // Update norms
-        if (k < nIter - 1)
+        if (k < sz - 1)
         {
             label q = k + 1;
-            for (const auto& val : RMatrix(A.subRow(k, q)))
+            for (const auto& val : RMatrix(AT.subColumn(k, q)))
             {
-                colNorms[q] -= magSqr(val);
+                norms[q] -= magSqr(val);
                 ++q;
             }
         }
     }
 
-    if (outputType_ == outputTypes::REDUCED_R)
+    calcQ(AT);
+
+    calcR(AT, Rdiag);
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::calcQ
+(
+    const MatrixType& AT
+)
+{
+    if (output_ == outputs::ONLY_R)
     {
-        A.resize(A.n(), A.n());
+        return;
+    }
+
+    const label n = AT.m();
+    const label m = AT.n();
+
+    Q_.resize(m, sz_);
+
+    MatrixType QT(Q_.transpose());
+
+    for (label k = sz_ - 1; k >= 0; --k)
+    {
+        for (label i = 0; i < m; ++i)
+        {
+            QT(k,i) = Zero;
+        }
+
+        QT(k,k) = scalar(1);
+
+        for (label j = k; j < sz_; ++j)
+        {
+            if (n > k && mag(AT(k,k)) > SMALL)
+            {
+                cmptType s(Zero);
+
+                for (label i = k; i < m; ++i)
+                {
+                    s += AT(k,i)*QT(j,i);
+                }
+
+                s /= -AT(k,k);
+
+                for (label i = k; i < m; ++i)
+                {
+                    QT(j,i) += s*Detail::conj(AT(k,i));
+                }
+            }
+        }
+    }
+
+    Q_ = QT.T();
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::calcR
+(
+    const MatrixType& AT,
+    const List<cmptType>& diag
+)
+{
+    if (output_ == outputs::ONLY_Q)
+    {
+        return;
+    }
+
+    const label n = AT.m();
+
+    R_.resize(sz_, n);
+
+    for (label i = 0; i < sz_; ++i)
+    {
+        for (label j = 0; j < n; ++j)
+        {
+            if (i < j)
+            {
+                R_(i,j) = AT(j,i);
+            }
+            else if (i == j)
+            {
+                R_(i,j) = diag[i];
+            }
+        }
     }
 }
 
@@ -182,49 +378,36 @@ 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
+    const modes mode,
+    const outputs output,
+    const bool pivoting,
+    MatrixType& A
 )
 :
-    outputType_(outputType),
-    storeMethod_(storeMethod),
-    colPivot_(colPivot),
+    mode_(mode),
+    output_(output),
+    sz_(calcShapeFactor(A)),
     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_()
+    p_()
 {
-    decompose(A);
+    // Non-conjugate transpose of input matrix
+    MatrixType AT(A.transpose());
+    A.clear();
+
+    if (pivoting)
+    {
+        decompose(AT, pivoting);
+    }
+    else
+    {
+        decompose(AT);
+    }
+
+    AT.clear();
 }
 
 
@@ -232,127 +415,35 @@ template<class MatrixType>
 Foam::QRMatrix<MatrixType>::QRMatrix
 (
     const MatrixType& A,
-    const outputTypes outputType,
-    const colPivoting colPivot
+    const modes mode,
+    const outputs output,
+    const bool pivoting
 )
 :
-    outputType_(outputType),
-    storeMethod_(storeMethods::OUT_OF_PLACE),
-    colPivot_(colPivot),
+    mode_(mode),
+    output_(output),
+    sz_(calcShapeFactor(A)),
     Q_(),
     R_(),
-    P_()
+    p_()
 {
-    decompose(A);
+    MatrixType AT(A.transpose());
+
+    if (pivoting)
+    {
+        decompose(AT, pivoting);
+    }
+    else
+    {
+        decompose(AT);
+    }
+
+    AT.clear();
 }
 
 
 // * * * * * * * * * * * * * * 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
 (
@@ -407,6 +498,62 @@ Foam::QRMatrix<MatrixType>::solve
 }
 
 
+template<class MatrixType>
+Foam::RectangularMatrix<typename MatrixType::cmptType>
+Foam::QRMatrix<MatrixType>::solve
+(
+    const RMatrix& rhs
+)
+{
+    const label mRows = R_.m();
+    const label nCols = R_.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(R_.diag());
+        for (const auto& val : diag)
+        {
+            if (mag(val) < SMALL)
+            {
+                WarningInFunction
+                    << "SMALL-valued diagonal elem in back-substitution."
+                    << endl;
+            }
+        }
+    }
+    #endif
+
+    RMatrix b({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 -= b(k, i)*R_(j, k);
+            }
+
+            b(j, i) = alpha/R_(j, j);
+        }
+    }
+
+    return b;
+}
+
+
 template<class MatrixType>
 typename Foam::QRMatrix<MatrixType>::SMatrix
 Foam::QRMatrix<MatrixType>::inv() const
@@ -430,72 +577,15 @@ 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
+MatrixType Foam::MatrixTools::pinv
 (
     const MatrixType& A,
-    const scalar tolerance
+    scalar tol
 )
 {
-    scalar tol = tolerance;
     typedef typename MatrixType::cmptType cmptType;
 
     if (A.empty())
@@ -518,15 +608,16 @@ MatrixType Foam::pinv
     QRMatrix<MatrixType> QRM
     (
         A,
-        QRMatrix<MatrixType>::outputTypes::FULL_QR,
-        QRMatrix<MatrixType>::colPivoting::TRUE
+        QRMatrix<MatrixType>::modes::FULL,
+        QRMatrix<MatrixType>::outputs::BOTH_QR,
+        QRMatrix<MatrixType>::pivoting::TRUE
     );
-    const MatrixType& R(QRM.R());
-    const MatrixType& Q(QRM.Q());
+    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;
+    label elemi = 0;
     label i = 0;
     while (i < 2)
     {
@@ -534,14 +625,14 @@ MatrixType Foam::pinv
 
         auto lessThan = [=](const cmptType& x) { return tol > mag(x); };
 
-        firstZeroElemi =
+        elemi =
             std::distance
             (
                 diag.cbegin(),
                 std::find_if(diag.cbegin(), diag.cend(), lessThan)
             );
 
-        if (firstZeroElemi == 0)
+        if (elemi == 0)
         {
             if (i == 0)
             {
@@ -572,7 +663,7 @@ MatrixType Foam::pinv
 
     // 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));
+    const RectangularMatrix<cmptType> R1(R.subMatrix(0, 0, elemi));
 
     // R2 (KP:p. 648)
     if (R1.n() < R1.m())
@@ -582,15 +673,15 @@ MatrixType Foam::pinv
         QRMatrix<SquareMatrix<cmptType>> QRSolve
         (
             C,
-            QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+            QRMatrix<SquareMatrix<cmptType>>::modes::FULL,
+            QRMatrix<SquareMatrix<cmptType>>::outputs::BOTH_QR
         );
 
         RectangularMatrix<cmptType> R2
         (
-            QRSolve.backSubstitution
+            QRSolve.solve
             (
-                QRSolve.R(),
-                QRSolve.Q() & (R1.T())
+                QRSolve.Q() & R1.T()
             )
         );
 
@@ -606,14 +697,14 @@ MatrixType Foam::pinv
         QRMatrix<SquareMatrix<cmptType>> QRSolve
         (
             C,
-            QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+            QRMatrix<SquareMatrix<cmptType>>::modes::FULL,
+            QRMatrix<SquareMatrix<cmptType>>::outputs::BOTH_QR
         );
 
         RectangularMatrix<cmptType> R2
         (
-            QRSolve.backSubstitution
+            QRSolve.solve
             (
-                QRSolve.R(),
                 QRSolve.Q() & R1
             )
         );
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index ee9bd62edf..dcb32cc6db 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,29 +28,36 @@ Class
     Foam::QRMatrix
 
 Description
-    QRMatrix (i.e. QR decomposition, QR factorisation or orthogonal-triangular
-    decomposition) decomposes a scalar/complex matrix \c A into the following
-    matrix product:
+    \c QRMatrix computes QR decomposition of a given
+    scalar/complex matrix \c A into the following:
 
     \verbatim
-        A = Q*R,
+        A = Q R
+    \endverbatim
+
+    or in case of QR decomposition with column pivoting:
+
+    \verbatim
+        A P = Q R
     \endverbatim
 
     where
-     \c Q is a unitary similarity matrix,
-     \c R is an upper triangular matrix.
+    \vartable
+        Q | Unitary/orthogonal matrix
+        R | Upper triangular matrix
+        P | Permutation matrix
+    \endvartable
 
-    Reference:
+    References:
     \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
+        TNT implementation:
+            Pozo, R. (1997).
+            Template Numerical Toolkit for linear algebra:
+            High performance programming with C++
+            and the Standard Template Library.
+            The International Journal of Supercomputer Applications
+            and High Performance Computing, 11(3), 251-263.
+            DOI:10.1177/109434209701100307
 
         QR decomposition with column pivoting (tag:QSB):
             Quintana-Ortí, G., Sun, X., & Bischof, C. H. (1998).
@@ -75,54 +82,47 @@ Description
             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
+    \heading Input:
+    \vartable
+        A | RectangularMatrix<Type> or SquareMatrix<Type>
+    \endvartable
 
-    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)
+    \heading Options for the decomposition mode:
+    \vartable
+        modes::FULL     | compute full-size decomposition
+        modes::ECONOMY  | compute economy-size decomposition
+    \endvartable
 
-    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
+    \heading Options for the output types:
+    \vartable
+        outputs::ONLY\_Q     | compute only Q
+        outputs::ONLY\_R     | compute only R
+        outputs::BOTH\_QR    | compute both Q and R
+    \endvartable
 
-    Options for the computation of column pivoting:
-     - colPivoting::FALSE:            switches off column pivoting
-     - colPivoting::TRUE:             switches on column pivoting
+    \heading Options for the column pivoting:
+    \vartable
+        pivoting::FALSE     | switch off column pivoting
+        pivoting::TRUE      | switch on column pivoting
+    \endvartable
 
-    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
+    \heading Output:
+    \vartable
+        Q | m-by-m (FULL) or m-by-k (ECONOMY) with k = min(m,n)
+        R | m-by-n (FULL) or k-by-n (ECONOMY) with k = min(m,n)
+        p | n-element label list
+        P | n-by-n permutation matrix
+    \endvartable
 
 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.
+    - \c QRMatrix involves modified implementations of the public-domain
+    library \c TNT, which is available at https://math.nist.gov/tnt/index.html.
+    - \c Q and \c R are always of the same type of the matrix \c A.
+    - \c Type can be \c scalar or \c complex.
 
 See also
     Test-QRMatrix.C
@@ -133,12 +133,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef QRMatrix_H
-#define QRMatrix_H
+#ifndef Foam_QRMatrix_H
+#define Foam_QRMatrix_H
 
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "complex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -146,39 +145,38 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                            Class QRMatrix Declaration
+                          Class QRMatrix Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class MatrixType>
 class QRMatrix
 {
-
 public:
 
     typedef typename MatrixType::cmptType cmptType;
     typedef SquareMatrix<cmptType> SMatrix;
     typedef RectangularMatrix<cmptType> RMatrix;
 
-    //- Options for the output matrix forms of QRMatrix
-    enum outputTypes : uint8_t
+    //- Options for the decomposition mode
+    enum modes : 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
+        FULL = 1,         //!< compute full-size decomposition
+        ECONOMY = 2,      //!< compute economy-size decomposition
     };
 
-    //- Options where to store R
-    enum storeMethods : uint8_t
+    //- Options for the output types
+    enum outputs : uint8_t
     {
-        IN_PLACE = 1,       //!< replaces input matrix content with \c R
-        OUT_OF_PLACE = 2    //!< creates new object of \c R
+        ONLY_Q = 1,      //!< compute only Q
+        ONLY_R = 2,      //!< compute only R
+        BOTH_QR = 3      //!< compute both Q and R
     };
 
-    //- Options for the computation of column pivoting
-    enum colPivoting : bool
+    //- Options for the column pivoting
+    enum pivoting : bool
     {
-        FALSE = false,      //!< switches off column pivoting
-        TRUE  = true        //!< switches on column pivoting
+        FALSE = false,      //!< switch off column pivoting
+        TRUE  = true        //!< switch on column pivoting
     };
 
 
@@ -186,53 +184,50 @@ private:
 
     // Private Data
 
-        //- Selected option for the output matrix forms of QRMatrix
-        outputTypes outputType_;
+        //- Decomposition mode
+        const modes mode_;
 
-        //- Selected option where to store R
-        const storeMethods storeMethod_;
+        //- Output type
+        const outputs output_;
 
-        //- Selected option for the computation of column pivoting
-        const colPivoting colPivot_;
+        //- Output shape factor based on the decomposition mode
+        const label sz_;
 
-        //- Unitary similarity matrix
+        //- Unitary/orthogonal matrix
         MatrixType Q_;
 
         //- Upper triangular matrix
         MatrixType R_;
 
         //- Permutation list (if column-pivoting is on)
-        labelList P_;
+        labelList p_;
 
 
     // Private Member Functions
 
-        //- Compute the QR decomposition without the column pivoting
-        void qr
-        (
-            MatrixType& A
-        );
+    // Evaluation
 
-        //- Compute the QR decomposition with the column pivoting
-        //  (QSB:Fig. 1)
-        void qrPivot
-        (
-            MatrixType& A
-        );
+        //- Calculate output shape factor based on the decomposition mode
+        label calcShapeFactor(const MatrixType& A) const;
 
-        //- Compute output matrices in selected forms
-        //- using Householder reflectors
-        //  (M:Section 4)
-        inline void applyHouseholder
-        (
-            MatrixType& A,
-            const RMatrix& reflector,
-            const label k = 0
-        );
+        //- Compute QR decomposition
+        void decompose(MatrixType& AT);
 
-        //- 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
+        //- Compute QR decomposition with column pivoting
+        void decompose(MatrixType& AT, const bool pivot);
+
+        //- Calculate Q
+        void calcQ(const MatrixType& AT);
+
+        //- Calculate R
+        void calcR(const MatrixType& AT, const List<cmptType>& diag);
+
+
+    // Linear system solution
+
+        //- 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
         template<template<typename> class ListContainer>
         void solvex(ListContainer<cmptType>& x) const;
 
@@ -246,33 +241,11 @@ private:
         ) const;
 
 
-protected:
+        //- No copy construct
+        QRMatrix(const QRMatrix&) = delete;
 
-    // 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
-        );
+        //- No copy assignment
+        QRMatrix& operator=(const QRMatrix&) = delete;
 
 
 public:
@@ -280,65 +253,66 @@ public:
     // Constructors
 
         //- Construct null
-        QRMatrix();
+        QRMatrix() = delete;
 
-        //- Construct QRMatrix without performing the decomposition
-        QRMatrix
+        //- Construct with a matrix and perform QR decomposition
+        explicit QRMatrix
         (
-            const outputTypes outputType,
-            const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
-            const colPivoting colPivot = colPivoting::FALSE
+            const modes mode,
+            const outputs output,
+            const bool pivoting,
+            MatrixType& A
         );
 
-        //- 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
+        //- Construct with a const matrix and perform QR decomposition
+        explicit QRMatrix
         (
             const MatrixType& A,
-            const outputTypes outputType,
-            const colPivoting colPivot = colPivoting::FALSE
+            const modes mode,
+            const outputs output = outputs::BOTH_QR,
+            const bool pivoting = false
         );
 
 
     // Member Functions
 
-        // Information
+    // Access
 
-        //- Return the unitary similarity matrix
-        //  Includes implicit round-to-zero as mutable operation
-        inline const MatrixType& Q() const;
+        //- Return const reference to Q
+        const MatrixType& Q() const noexcept
+        {
+            return Q_;
+        }
 
-        //- Return the upper triangular matrix
-        inline const MatrixType& R() const;
+        //- Return reference to Q
+        MatrixType& Q() noexcept
+        {
+            return Q_;
+        }
 
-        //- Return the permutation order (P) as a list
-        inline const labelList& orderP() const;
+        //- Return const reference to R
+        const MatrixType& R() const noexcept
+        {
+            return R_;
+        }
+
+        //- Return reference to R
+        MatrixType& R() noexcept
+        {
+            return R_;
+        }
+
+        //- Return const reference to p
+        const labelList& p() const noexcept
+        {
+            return p_;
+        }
 
         //- 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
-        );
+    // Linear system solution
 
         //- Solve the linear system with the given source
         //- and return the solution in the argument x
@@ -372,41 +346,39 @@ public:
             const IndirectListBase<cmptType, Addr>& source
         ) const;
 
-        //- Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source
-        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
+        //- Solve a row-echelon-form linear system (Ax = b)
+        //- starting from the bottom by back substitution
+        //  A = R: Non-singular upper-triangular square matrix (m-by-m)
+        //  b: Source (m-by-p)
+        //  x: Solution (m-by-p)
+        RMatrix solve
         (
-            const SMatrix& A,
             const RMatrix& b
         );
+
+        //- Return the inverse of (Q*R), solving x = (Q*R).inv()*source
+        SMatrix inv() const;
 };
 
+namespace MatrixTools
+{
 
 // * * * * * * * * * * * * * * 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.
+//- Moore-Penrose inverse of singular/non-singular square/rectangular
+//- scalar/complex matrices (KPP:p. 9834; KP:p. 648)
+//  The tolerance to ensure the R1 matrix full-rank is set to 1e-5
+//  by (TA; mentioned in (KPP:p. 9832)) in contrast to 1e-13 (KPP:p. 9834).
 template<class MatrixType>
 MatrixType pinv
 (
     const MatrixType& A,
-    const scalar tol = 1e-5
+    scalar tol = 1e-5
 );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace MatrixTools
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
index 929f8b6fc8..bf61a6b4e9 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,166 +26,20 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class MatrixType>
-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 MatrixType& Foam::QRMatrix<MatrixType>::Q() const
-{
-    const_cast<MatrixType&>(Q_).round();
-
-    return Q_;
-}
-
-
-template<class MatrixType>
-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>
+inline Foam::SquareMatrix<typename MatrixType::cmptType>
 Foam::QRMatrix<MatrixType>::P() const
 {
-    SquareMatrix<cmptType> permMat(P_.size(), Zero);
+    SquareMatrix<cmptType> P(p_.size(), Zero);
 
-    forAll(P_, jcol)
+    forAll(p_, jcol)
     {
-        permMat(P_[jcol], jcol) = pTraits<cmptType>::one;
+        P(p_[jcol], jcol) = pTraits<cmptType>::one;
     }
 
-    return permMat;
+    return P;
 }