LLTMatrix: New matrix form to support Cholesky decomposition

of symmetric positive-definite matrices and the solution of associated
linear systems.

applications/test/Matrix/Test-Matrix.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "scalarMatrices.H"
+#include "LLTMatrix.H"
 #include "vector.H"
 #include "tensor.H"
 #include "IFstream.H"
@@ -139,6 +140,27 @@ int main(int argc, char *argv[])
         Info<< "det = " << detDecomposed(squareMatrix, sign) << 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);
+
+        LLTMatrix<scalar> LLT(squareMatrix);
+        scalarField x(LLT.solve(source));
+
+        Info<< "LLT solve residual " << (squareMatrix*x - source) << endl;
+    }
+
     Info<< "\nEnd\n" << endl;
 
     return 0;

src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LLTMatrix.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::LLTMatrix<Type>::decompose()
+{
+    SquareMatrix<Type>& LLT = *this;
+    const label m = LLT.m();
+
+    for (label i = 0; i < m; i++)
+    {
+        for (label j = 0; j < m; j++)
+        {
+            if (j > i)
+            {
+                LLT(i, j) = Zero;
+                continue;
+            }
+
+            Type sum = LLT(i, j);
+
+            for (label k = 0; k < j; k++)
+            {
+                sum -= LLT(i, k)*LLT(j, k);
+            }
+
+            if (i > j)
+            {
+                LLT(i, j) = sum/LLT(j, j);
+            }
+            else if (sum > 0)
+            {
+                LLT(i, i) = sqrt(sum);
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Cholesky decomposition failed, "
+                       "matrix is not symmetric positive definite"
+                    << abort(FatalError);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::LLTMatrix<Type>::LLTMatrix(const SquareMatrix<Type>& M)
+:
+    SquareMatrix<Type>(M)
+{
+    decompose();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::LLTMatrix<Type>::solve
+(
+    Field<Type>& x,
+    const Field<Type>& source
+) const
+{
+    const SquareMatrix<Type>& LLT = *this;
+    const label m = LLT.m();
+
+    for (label i = 0; i < m; i++)
+    {
+        Type sum = source[i];
+
+        for (label j = 0; j < i; j++)
+        {
+            sum = sum - LLT(i, j)*x[j];
+        }
+
+        x[i] = sum/LLT(i, i);
+    }
+
+    for (int i = m - 1; i >= 0; i--)
+    {
+        Type sum = x[i];
+
+        for (label j = i + 1; j < m; j++)
+        {
+            sum = sum - LLT(j, i)*x[j];
+        }
+
+        x[i] = sum/LLT(i, i);
+    }
+
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::LLTMatrix<Type>::solve
+(
+    const Field<Type>& source
+) const
+{
+    tmp<Field<Type>> tx(new Field<Type>(this->m()));
+    Field<Type>& x = tx.ref();
+
+    solve(x, source);
+
+    return tx;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::LLTMatrix
+
+Description
+    Templated class to perform the Cholesky decomposition on a
+    symmetric positive-definite matrix.
+
+    Member functions are provided to solve linear systems using the LLT
+    decomposition.
+
+SourceFiles
+    LLTMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LLTMatrix_H
+#define LLTMatrix_H
+
+#include "SquareMatrix.H"
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class LLTMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class LLTMatrix
+:
+    public SquareMatrix<Type>
+{
+    // Private Member Functions
+
+        //- Perform the Cholesky decomposition of the matrix
+        void decompose();
+
+public:
+
+    // Constructors
+
+        //- Construct from a square matrix and perform the decomposition
+        LLTMatrix(const SquareMatrix<Type>& M);
+
+
+    // Member Functions
+
+        //- Solve the linear system with the given source
+        //  and returning the solution in the Field argument x
+        void solve(Field<Type>& x, const Field<Type>& source) const;
+
+        //- Solve the linear system with the given source
+        //  returning the solution
+        tmp<Field<Type>> solve(const Field<Type>& source) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "LLTMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //