From e9199c6e14ddef7b19d0d791b3986701ecba86f9 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 24 Mar 2016 14:49:25 +0000
Subject: [PATCH] QRMatrix: New class to provide QR-decomposition by
 Householder reflection

This development is sponsored by Carnegie Wave Energy Ltd.
---
 applications/test/Matrix/Test-Matrix.C     |  27 +++
 src/OpenFOAM/matrices/QRMatrix/QRMatrix.C  | 249 +++++++++++++++++++++
 src/OpenFOAM/matrices/QRMatrix/QRMatrix.H  | 133 +++++++++++
 src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H |  51 +++++
 4 files changed, 460 insertions(+)
 create mode 100644 src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
 create mode 100644 src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
 create mode 100644 src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H

diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index c8691e81b5..87a39386fc 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -25,6 +25,7 @@ License
 
 #include "scalarMatrices.H"
 #include "LLTMatrix.H"
+#include "QRMatrix.H"
 #include "vector.H"
 #include "tensor.H"
 #include "IFstream.H"
@@ -161,6 +162,32 @@ int main(int argc, char *argv[])
         Info<< "LLT solve residual " << (squareMatrix*x - source) << endl;
     }
 
+    {
+        scalarSquareMatrix squareMatrix(3, Zero);
+
+        squareMatrix(0, 0) = 4;
+        squareMatrix(0, 1) = 12;
+        squareMatrix(0, 2) = -16;
+        squareMatrix(1, 0) = 12;
+        squareMatrix(1, 1) = 37;
+        squareMatrix(1, 2) = -43;
+        squareMatrix(2, 0) = -16;
+        squareMatrix(2, 1) = -43;
+        squareMatrix(2, 2) = 98;
+
+        scalarField source(3, 1);
+
+        QRMatrix<scalarSquareMatrix> QR(squareMatrix);
+        scalarField x(QR.solve(source));
+
+        Info<< "QR solve residual "
+            << (squareMatrix*x - source) << endl;
+
+        Info<< "QR inverse solve residual "
+            << (x - QR.inverse()*source) << endl;
+
+    }
+
     Info<< "\nEnd\n" << endl;
 
     return 0;
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
new file mode 100644
index 0000000000..80e33b0fb2
--- /dev/null
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -0,0 +1,249 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "QRMatrix.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix(const MatrixType& M)
+{
+    decompose(M);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose(const MatrixType& M)
+{
+    const label m = M.m();
+    const label n = M.n();
+
+    // Initialize the R-matrix to M
+    R_ = M;
+
+    // Initialize the Q-matrix to I
+    Q_.setSize(m);
+    Q_ = I;
+
+    // Pre-allocate temporary storage for the Householder steps
+    QMatrixType Qk(m);
+    QMatrixType Rk(m);
+    Field<cmptType> uk(m);
+
+    for (label k=0; k<n; k++)
+    {
+        // alpha = -|column k of Rk|
+        cmptType alpha = Zero;
+        for (label bi=k; bi<m; bi++)
+        {
+            alpha += sqr(R_(bi, k));
+        }
+        alpha = sqrt(alpha);
+
+        if (R_(k, k) > 0)
+        {
+            alpha = -alpha;
+        }
+
+        // uk = column k of Rk - alpha*ek
+        // rSumSqrUk = 2/sum(sqr(uk))
+        uk[k] = R_(k, k) - alpha;
+        cmptType rSumSqrUk = sqr(uk[k]);
+        for (label bi=k+1; bi<m; bi++)
+        {
+            uk[bi] = R_(bi, k);
+            rSumSqrUk += sqr(uk[bi]);
+        }
+        rSumSqrUk = 2/rSumSqrUk;
+
+        // Qk = I - 2*u*uT/sum(sqr(uk))
+        for (label bi=k; bi<m; bi++)
+        {
+            for (label bj=k; bj<m; bj++)
+            {
+                Qk(bi, bj) = -rSumSqrUk*uk[bi]*uk[bj];
+            }
+            Qk(bi, bi) += 1;
+        }
+
+        // Rk = Qk*R
+        for (label bi=k; bi<m; bi++)
+        {
+            for (label bk=k; bk<n; bk++)
+            {
+                Rk(bi, bk) = Zero;
+                for (label bj=k; bj<n; bj++)
+                {
+                    Rk(bi, bk) += Qk(bi, bj)*R_(bj, bk);
+                }
+            }
+        }
+
+        // Ensure diagonal is positive
+        if (R_(k, k) < 0)
+        {
+            // R = -Rk
+            // Q = -Q
+            for (label bi=k; bi<m; bi++)
+            {
+                for (label bj=k; bj<n; bj++)
+                {
+                    R_(bi, bj) = -Rk(bi, bj);
+                }
+                for (label bj=k; bj<m; bj++)
+                {
+                    Q_(bi, bj) = -Q_(bi, bj);
+                }
+            }
+        }
+        else
+        {
+            // R = Rk
+            for (label bi=k; bi<m; bi++)
+            {
+                for (label bj=k; bj<n; bj++)
+                {
+                    R_(bi, bj) = Rk(bi, bj);
+                }
+            }
+        }
+
+        // Q = Q*Qk (using Rk as temporary storage)
+        for (label bi=0; bi<m; bi++)
+        {
+            for (label bk=k; bk<m; bk++)
+            {
+                Rk(bi, bk) = Zero;
+                for (label bj=k; bj<m; bj++)
+                {
+                    Rk(bi, bk) += Q_(bi, bj)*Qk(bj, bk);
+                }
+            }
+        }
+        for (label bi=0; bi<m; bi++)
+        {
+            for (label bj=k; bj<n; bj++)
+            {
+                Q_(bi, bj) = Rk(bi, bj);
+            }
+        }
+    }
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::solvex
+(
+    Field<cmptType>& x
+) const
+{
+    const label n = R_.n();
+
+    for (int i=n-1; i>=0; i--)
+    {
+        cmptType sum = x[i];
+
+        for (label j=i+1; j<n; j++)
+        {
+            sum -= x[j]*R_(i, j);
+        }
+
+        if (mag(R_(i, i)) < SMALL)
+        {
+            FatalErrorInFunction
+                << "Back-substitution failed due to small diagonal"
+                << abort(FatalError);
+        }
+
+        x[i] = sum/R_(i, i);
+    }
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::solve
+(
+    Field<cmptType>& x,
+    const Field<cmptType>& source
+) const
+{
+    const label m = Q_.m();
+
+    // x = Q_.T()*source;
+    for (label i=0; i<m; i++)
+    {
+        x[i] = 0;
+        for (label j=0; j<m; j++)
+        {
+            x[i] += Q_(j, i)*source[j];
+        }
+    }
+
+    solvex(x);
+}
+
+
+template<class MatrixType>
+Foam::tmp<Foam::Field<typename MatrixType::cmptType>>
+Foam::QRMatrix<MatrixType>::solve
+(
+    const Field<cmptType>& source
+) const
+{
+    tmp<Field<cmptType>> tx(new Field<cmptType>(Q_.m()));
+    Field<cmptType>& x = tx.ref();
+
+    solve(x, source);
+
+    return tx;
+}
+
+
+template<class MatrixType>
+typename Foam::QRMatrix<MatrixType>::QMatrixType
+Foam::QRMatrix<MatrixType>::inverse() const
+{
+    const label m = Q_.m();
+
+    Field<cmptType> x(m);
+    QMatrixType inv(m);
+
+    for (label i=0; i<m; i++)
+    {
+        for (label j=0; j<m; j++)
+        {
+            x[j] = Q_(i, j);
+        }
+        solvex(x);
+        inv.block(m, 1, 0, i) = x;
+    }
+
+    return inv;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
new file mode 100644
index 0000000000..bbad659be2
--- /dev/null
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::QRMatrix
+
+Description
+    Class templated on matrix type to perform the QR decomposition using
+    Householder reflections on a square or rectangular matrix.
+
+SourceFiles
+    QRMatrixI.H
+    QRMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef QRMatrix_H
+#define QRMatrix_H
+
+#include "SquareMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class QRMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MatrixType>
+class QRMatrix
+{
+
+public:
+
+    typedef typename MatrixType::cmptType cmptType;
+    typedef SquareMatrix<cmptType> QMatrixType;
+    typedef MatrixType RMatrixType;
+
+private:
+
+    // Private data
+
+        //- The Q-matrix
+        QMatrixType Q_;
+
+        //- The R-matrix
+        RMatrixType R_;
+
+
+    // Private member functions
+
+        //- 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
+        void solvex(Field<cmptType>& x) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline QRMatrix();
+
+        //- Construct decomposing given matrix
+        QRMatrix(const MatrixType& M);
+
+
+    // Member Functions
+
+        //- Return Q-matrix
+        inline const QMatrixType& Q() const;
+
+        //- Return R-matrix
+        inline const RMatrixType& R() const;
+
+        //- Decompose given matrix
+        void decompose(const MatrixType& M);
+
+        //- Solve the linear system with the given source
+        //  and returning the solution in the Field argument x
+        void solve(Field<cmptType>& x, const Field<cmptType>& source) const;
+
+        //- Solve the linear system with the given source
+        //  returning the solution
+        tmp<Field<cmptType>> solve(const Field<cmptType>& source) const;
+
+        //- Return the inverse of a square matrix
+        QMatrixType inverse() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "QRMatrixI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "QRMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
new file mode 100644
index 0000000000..1d9ca4b620
--- /dev/null
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+inline Foam::QRMatrix<MatrixType>::QRMatrix()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MatrixType>
+inline const typename Foam::QRMatrix<MatrixType>::QMatrixType&
+Foam::QRMatrix<MatrixType>::Q() const
+{
+    return Q_;
+}
+
+
+template<class MatrixType>
+inline const typename Foam::QRMatrix<MatrixType>::RMatrixType&
+Foam::QRMatrix<MatrixType>::R() const
+{
+    return R_;
+}
+
+
+// ************************************************************************* //