LLTMatrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "LLTMatrix.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class Type>
32 {}
33 
34 
35 template<class Type>
37 {
38  decompose(M);
39 }
40 
41 
42 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
43 
44 template<class Type>
46 {
47  SquareMatrix<Type>& LLT = *this;
48 
49  // Initialise the LLT decomposition matrix to M
50  LLT = M;
51 
52  const label m = LLT.m();
53 
54  for (label i=0; i<m; i++)
55  {
56  for (label j=0; j<m; j++)
57  {
58  if (j > i)
59  {
60  LLT(i, j) = Zero;
61  continue;
62  }
63 
64  Type sum = LLT(i, j);
65 
66  for (label k=0; k<j; k++)
67  {
68  sum -= LLT(i, k)*LLT(j, k);
69  }
70 
71  if (i > j)
72  {
73  LLT(i, j) = sum/LLT(j, j);
74  }
75  else if (sum > 0)
76  {
77  LLT(i, i) = sqrt(sum);
78  }
79  else
80  {
82  << "Cholesky decomposition failed, "
83  "matrix is not symmetric positive definite"
84  << abort(FatalError);
85  }
86  }
87  }
88 }
89 
90 
91 template<class Type>
93 (
94  Field<Type>& x,
95  const Field<Type>& source
96 ) const
97 {
98  // If x and source are different initialise x = source
99  if (&x != &source)
100  {
101  x = source;
102  }
103 
104  const SquareMatrix<Type>& LLT = *this;
105  const label m = LLT.m();
106 
107  for (label i=0; i<m; i++)
108  {
109  Type sum = source[i];
110 
111  for (label j=0; j<i; j++)
112  {
113  sum = sum - LLT(i, j)*x[j];
114  }
115 
116  x[i] = sum/LLT(i, i);
117  }
118 
119  for (int i=m - 1; i >= 0; i--)
120  {
121  Type sum = x[i];
122 
123  for (label j=i + 1; j<m; j++)
124  {
125  sum = sum - LLT(j, i)*x[j];
126  }
127 
128  x[i] = sum/LLT(i, i);
129  }
130 
131 }
132 
133 
134 template<class Type>
136 (
137  const Field<Type>& source
138 ) const
139 {
140  tmp<Field<Type>> tx(new Field<Type>(this->m()));
141  Field<Type>& x = tx.ref();
142 
143  solve(x, source);
144 
145  return tx;
146 }
147 
148 
149 // ************************************************************************* //
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedScalar sqrt(const dimensionedScalar &ds)
label k
Boltzmann constant.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Pre-declare SubField and related Field type.
Definition: Field.H:56
void solve(Field< Type > &x, const Field< Type > &source) const
Solve the linear system with the given source.
Definition: LLTMatrix.C:93
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
LLTMatrix()
Construct null.
Definition: LLTMatrix.C:31
rhoEqn solve()
label m() const
Return the number of rows.
Definition: MatrixI.H:57
A class for managing temporary objects.
Definition: PtrList.H:53
void decompose(const SquareMatrix< Type > &M)
Perform the Cholesky decomposition of the matrix.
Definition: LLTMatrix.C:45
#define M(I)
A templated 2D square matrix of objects of <T>, where the n x n matrix dimension is known and used fo...
Definition: SquareMatrix.H:58