ENH: QRMatrix: refactor the QR decomposition algorithms (fixes #2212)

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;
+        }
+    }
+}
+
+
 // ************************************************************************* //

applications/test/matrices/QRMatrix/Make/options

@@ -1,2 +1,2 @@
-/* EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude */
+EXE_INC = -I../../TestTools
 /* EXE_LIBS = -lfiniteVolume */

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());
+    if (mode_ == modes::FULL)
-
-    // Loop through all subcolumns of which the diagonal elem is the first elem
-    for (label k = 0; k < nIter; ++k)
     {
-        const RMatrix u(householderReflector(A.subColumn(k, k)));
+        return A.m();
-
+    }
-        applyHouseholder(A, u, k);
+    else if (mode_ == modes::ECONOMY)
+    {
+        return min(A.m(), A.n());
     }
 
-    if (outputType_ == outputTypes::REDUCED_R)
+    return 0;
-    {
-        A.resize(A.n(), A.n());
-    }
 }
 
 
 template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::qrPivot
+void Foam::QRMatrix<MatrixType>::decompose
 (
-    MatrixType& A
+    MatrixType& AT
 )
 {
-    const label nCols = A.n();
+    const label n = AT.m();
-    const label nIter = min(A.m() - 1, A.n());
+    const label m = AT.n();
 
-    // Initialise permutation vector, and column norm vector
+    List<cmptType> Rdiag(n, Zero);
-    P_ = identity(nCols);
 
-    // Initialise vector norms of each column of A
+    for (label k = 0; k < n; ++k)
-    List<scalar> colNorms(nCols);
-    for (label k = 0; k < nCols; ++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
+    calcQ(AT);
-    for (label k = 0; k < nIter; ++k)
-    {
-        const labelRange colRange(k, nCols);
-        const SubList<scalar> subColNorms(colNorms, colRange);
 
-        // Column pivoting
+    calcR(AT, Rdiag);
-        const label maxColNormi =
+}
+
+
+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(),
+                it,
-                std::max_element(subColNorms.cbegin(), subColNorms.cend())
+                std::max_element(it, norms.cend())
             );
 
-        // Swap R_, P_ and colNorms_ according to pivot column if the current
+        // Swap A, P and norms according to pivot column if the current
-        // column is not the max norm column by avoiding any column swap where
+        // column is not the max norm column. Also avoid any column swaps
-        // the leading elem is very small
+        // where the leading element is very small
-        if (maxColNormi != 0 && SMALL < mag(A(k, k + maxColNormi)))
+        if (mag(AT(k + maxNormi,k)) > SMALL && maxNormi != 0)
         {
-            const RMatrix R1(A.subColumn(k));
+            const RMatrix R1(AT.subRow(k));
-            const RMatrix R2(A.subColumn(maxColNormi + k));
+            const RMatrix R2(AT.subRow(maxNormi + k));
-            A.subColumn(k) = R2;
-            A.subColumn(maxColNormi + k) = R1;
 
-            Swap(P_[k], P_[maxColNormi + k]);
+            AT.subRow(k) = R2;
-            Swap(colNorms[k], colNorms[maxColNormi + k]);
+            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 modes mode,
-    const storeMethods storeMethod,
+    const outputs output,
-    const colPivoting colPivot
+    const bool pivoting,
+    MatrixType& A
 )
 :
-    outputType_(outputType),
+    mode_(mode),
-    storeMethod_(storeMethod),
+    output_(output),
-    colPivot_(colPivot),
+    sz_(calcShapeFactor(A)),
     Q_(),
     R_(),
-    P_()
+    p_()
-{}
-
-
-template<class MatrixType>
-Foam::QRMatrix<MatrixType>::QRMatrix
-(
-    MatrixType& A,
-    const outputTypes outputType,
-    const storeMethods storeMethod,
-    const colPivoting colPivot
-)
-:
-    outputType_(outputType),
-    storeMethod_(storeMethod),
-    colPivot_(colPivot),
-    Q_(),
-    R_(),
-    P_()
 {
-    decompose(A);
+    // 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 modes mode,
-    const colPivoting colPivot
+    const outputs output,
+    const bool pivoting
 )
 :
-    outputType_(outputType),
+    mode_(mode),
-    storeMethod_(storeMethods::OUT_OF_PLACE),
+    output_(output),
-    colPivot_(colPivot),
+    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>::modes::FULL,
-        QRMatrix<MatrixType>::colPivoting::TRUE
+        QRMatrix<MatrixType>::outputs::BOTH_QR,
+        QRMatrix<MatrixType>::pivoting::TRUE
     );
-    const MatrixType& R(QRM.R());
+    const MatrixType& R = QRM.R();
-    const MatrixType& Q(QRM.Q());
+    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
             )
         );

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
+    \c QRMatrix computes QR decomposition of a given
-    decomposition) decomposes a scalar/complex matrix \c A into the following
+    scalar/complex matrix \c A into the following:
-    matrix product:
 
     \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,
+    \vartable
-     \c R is an upper triangular matrix.
+        Q | Unitary/orthogonal matrix
+        R | Upper triangular matrix
+        P | Permutation matrix
+    \endvartable
 
-    Reference:
+    References:
     \verbatim
-        QR decomposition:
+        TNT implementation:
-            mathworld.wolfram.com/QRDecomposition.html (Retrieved:17-06-19)
+            Pozo, R. (1997).
-            w.wiki/52X (Retrieved: 17-06-19)
+            Template Numerical Toolkit for linear algebra:
-
+            High performance programming with C++
-        QR decomposition with Householder reflector (tag:M):
+            and the Standard Template Library.
-            Monahan, J. F. (2011).
+            The International Journal of Supercomputer Applications
-            Numerical methods of statistics.
+            and High Performance Computing, 11(3), 251-263.
-            Cambridge: Cambridge University Press.
+            DOI:10.1177/109434209701100307
-            DOI:10.1017/CBO9780511977176
 
         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:
+    \heading Input:
-     - \c Q is always of the type of the matrix \c A
+    \vartable
-     - \c R is always of the type of the matrix \c A
+        A | RectangularMatrix<Type> or SquareMatrix<Type>
+    \endvartable
 
-    Options for the output forms of \c QRMatrix (for an (m-by-n) input matrix
+    \heading Options for the decomposition mode:
-    \c A with k = min(m, n)):
+    \vartable
-     - outputTypes::FULL_R:     computes only \c R                   (m-by-n)
+        modes::FULL     | compute full-size decomposition
-     - outputTypes::FULL_QR:    computes both \c R and \c Q          (m-by-m)
+        modes::ECONOMY  | compute economy-size decomposition
-     - outputTypes::REDUCED_R:  computes only reduced \c R           (k-by-n)
+    \endvartable
 
-    Options where to store \c R:
+    \heading Options for the output types:
-     - storeMethods::IN_PLACE:        replaces input matrix content with \c R
+    \vartable
-     - storeMethods::OUT_OF_PLACE:    creates new object of \c R
+        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:
+    \heading Options for the column pivoting:
-     - colPivoting::FALSE:            switches off column pivoting
+    \vartable
-     - colPivoting::TRUE:             switches on column pivoting
+        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
+    \heading Output:
-    the following limitations:
+    \vartable
-     - \c A         = a scalar square matrix
+        Q | m-by-m (FULL) or m-by-k (ECONOMY) with k = min(m,n)
-     - output type  = outputTypes::FULL_QR
+        R | m-by-n (FULL) or k-by-n (ECONOMY) with k = min(m,n)
-     - store method = storeMethods::IN_PLACE
+        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.
+    - \c QRMatrix involves modified implementations of the public-domain
-    - The option combination:
+    library \c TNT, which is available at https://math.nist.gov/tnt/index.html.
-      - outputTypes::REDUCED_R
+    - \c Q and \c R are always of the same type of the matrix \c A.
-      - storeMethods::IN_PLACE
+    - \c Type can be \c scalar or \c complex.
-      will not modify the rows of input matrix \c A after its nth row.
-    - Both FULL_R and REDUCED_R QR decompositions execute the same number of
-      operations. Yet REDUCED_R QR decomposition returns only the first n rows
-      of \c R if m > n for an input m-by-n matrix \c A.
-    - For m <= n, FULL_R and REDUCED_R will produce the same matrices.
 
 See also
     Test-QRMatrix.C
@@ -133,12 +133,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef QRMatrix_H
+#ifndef Foam_QRMatrix_H
-#define 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
+    //- Options for the decomposition mode
-    enum outputTypes : uint8_t
+    enum modes : uint8_t
     {
-        FULL_R = 1,         //!< computes only \c R
+        FULL = 1,         //!< compute full-size decomposition
-        FULL_QR = 2,        //!< computes both \c R and \c Q
+        ECONOMY = 2,      //!< compute economy-size decomposition
-        REDUCED_R = 3       //!< computes only reduced \c R
     };
 
-    //- Options where to store R
+    //- Options for the output types
-    enum storeMethods : uint8_t
+    enum outputs : uint8_t
     {
-        IN_PLACE = 1,       //!< replaces input matrix content with \c R
+        ONLY_Q = 1,      //!< compute only Q
-        OUT_OF_PLACE = 2    //!< creates new object of \c R
+        ONLY_R = 2,      //!< compute only R
+        BOTH_QR = 3      //!< compute both Q and R
     };
 
-    //- Options for the computation of column pivoting
+    //- Options for the column pivoting
-    enum colPivoting : bool
+    enum pivoting : bool
     {
-        FALSE = false,      //!< switches off column pivoting
+        FALSE = false,      //!< switch off column pivoting
-        TRUE  = true        //!< switches on column pivoting
+        TRUE  = true        //!< switch on column pivoting
     };
 
 
@@ -186,53 +184,50 @@ private:
 
     // Private Data
 
-        //- Selected option for the output matrix forms of QRMatrix
+        //- Decomposition mode
-        outputTypes outputType_;
+        const modes mode_;
 
-        //- Selected option where to store R
+        //- Output type
-        const storeMethods storeMethod_;
+        const outputs output_;
 
-        //- Selected option for the computation of column pivoting
+        //- Output shape factor based on the decomposition mode
-        const colPivoting colPivot_;
+        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
+    // Evaluation
-        void qr
-        (
-            MatrixType& A
-        );
 
-        //- Compute the QR decomposition with the column pivoting
+        //- Calculate output shape factor based on the decomposition mode
-        //  (QSB:Fig. 1)
+        label calcShapeFactor(const MatrixType& A) const;
-        void qrPivot
-        (
-            MatrixType& A
-        );
 
-        //- Compute output matrices in selected forms
+        //- Compute QR decomposition
-        //- using Householder reflectors
+        void decompose(MatrixType& AT);
-        //  (M:Section 4)
-        inline void applyHouseholder
-        (
-            MatrixType& A,
-            const RMatrix& reflector,
-            const label k = 0
-        );
 
-        //- Solve the linear system with the Field argument x initialized to
+        //- Compute QR decomposition with column pivoting
-        //- the appropriate transformed source (e.g. Q.T()*source)
+        void decompose(MatrixType& AT, const bool pivot);
-        //- and return the solution in x
+
+        //- 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
+        //- No copy assignment
-
+        QRMatrix& operator=(const QRMatrix&) = delete;
-        //- Compute Householder reflector on a given matrix column, u
-        //  (M:Eq.6.814)
-        inline RMatrix householderReflector
-        (
-            RMatrix u
-        );
-
-        //- Apply (in-place) Householder reflectors from the left side: u*A
-        inline void applyLeftReflector
-        (
-            MatrixType& A,
-            const RMatrix& reflector,
-            const label k = 0,
-            const label k1 = 0
-        );
-
-        //- Apply (in-place) Householder reflectors from the right side: (u*A)*u
-        inline void applyRightReflector
-        (
-            MatrixType& A,
-            const RMatrix& reflector,
-            const label k = 0
-        );
 
 
 public:
@@ -280,65 +253,66 @@ public:
     // Constructors
 
         //- Construct null
-        QRMatrix();
+        QRMatrix() = delete;
 
-        //- Construct QRMatrix without performing the decomposition
+        //- Construct with a matrix and perform QR decomposition
-        QRMatrix
+        explicit QRMatrix
         (
-            const outputTypes outputType,
+            const modes mode,
-            const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
+            const outputs output,
-            const colPivoting colPivot = colPivoting::FALSE
+            const bool pivoting,
+            MatrixType& A
         );
 
-        //- Construct QRMatrix and perform the QR decomposition
+        //- Construct with a const matrix and perform QR decomposition
-        QRMatrix
+        explicit QRMatrix
-        (
-            MatrixType& A,
-            const outputTypes outputType,
-            const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
-            const colPivoting colPivot = colPivoting::FALSE
-        );
-
-        //- Construct QRMatrix and perform the QR decomposition
-        QRMatrix
         (
             const MatrixType& A,
-            const outputTypes outputType,
+            const modes mode,
-            const colPivoting colPivot = colPivoting::FALSE
+            const outputs output = outputs::BOTH_QR,
+            const bool pivoting = false
         );
 
 
     // Member Functions
 
-        // Information
+    // Access
 
-        //- Return the unitary similarity matrix
+        //- Return const reference to Q
-        //  Includes implicit round-to-zero as mutable operation
+        const MatrixType& Q() const noexcept
-        inline const MatrixType& Q() const;
+        {
+            return Q_;
+        }
 
-        //- Return the upper triangular matrix
+        //- Return reference to Q
-        inline const MatrixType& R() const;
+        MatrixType& Q() noexcept
+        {
+            return Q_;
+        }
 
-        //- Return the permutation order (P) as a list
+        //- Return const reference to R
-        inline const labelList& orderP() const;
+        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
+    // Linear system solution
-
-        //- Compute QR decomposition according to constructor settings
-        void decompose
-        (
-            MatrixType& A
-        );
-
-        //- Compute QR decomposition according to constructor settings
-        void decompose
-        (
-            const MatrixType& A
-        );
 
         //- Solve the linear system with the given source
         //- and return the solution in the argument x
@@ -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
+        //- Solve a row-echelon-form linear system (Ax = b)
-        SMatrix inv() const;
+        //- starting from the bottom by back substitution
-
+        //  A = R: Non-singular upper-triangular square matrix (m-by-m)
-        //- Solve a row-echelon-form linear system starting from the bottom
+        //  b: Source (m-by-p)
-        //  Back substitution: Ax = b where:
+        //  x: Solution (m-by-p)
-        //  A = Non-singular upper-triangular square matrix (m-by-m)
+        RMatrix solve
-        //  x = Solution (m-by-p)
-        //  b = Source (m-by-p)
-        RMatrix backSubstitution
         (
-            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
+//- Moore-Penrose inverse of singular/non-singular square/rectangular
-//- square/rectangular scalar/complex matrices
+//- scalar/complex matrices (KPP:p. 9834; KP:p. 648)
-//  (KPP:p. 9834; KP:p. 648)
+//  The tolerance to ensure the R1 matrix full-rank is set to 1e-5
-//  The tolerance (i.e. tol) to ensure the R1 matrix full-rank is given as 1e-13
+//  by (TA; mentioned in (KPP:p. 9832)) in contrast to 1e-13 (KPP:p. 9834).
-//  in the original paper (KPP:p. 9834).  Nonetheless, the tolerance = 1e-5
-//  given by (TA; mentioned in (KPP:p. 9832)) was observed to be more robust
-//  for the tested scenarios.
 template<class MatrixType>
 MatrixType pinv
 (
     const MatrixType& A,
-    const scalar tol = 1e-5
+    scalar tol = 1e-5
 );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace MatrixTools
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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
+inline Foam::SquareMatrix<typename MatrixType::cmptType>
-{
-    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>
 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;
 }