Nektar++
PreconditionerDiagonal.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: PreconditionerDiagonal.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: PreconditionerDiagonal definition
32//
33///////////////////////////////////////////////////////////////////////////////
34
39#include <cmath>
40
41using namespace std;
42
44{
45/**
46 * Registers the class with the Factory.
47 */
50 "Diagonal", PreconditionerDiagonal::create, "Diagonal Preconditioning");
51/**
52 * @class Preconditioner
53 *
54 * This class implements diagonal preconditioning for the conjugate
55 * gradient matrix solver.
56 */
57
59 const std::shared_ptr<GlobalLinSys> &plinsys,
60 const AssemblyMapSharedPtr &pLocToGloMap)
61 : Preconditioner(plinsys, pLocToGloMap)
62{
63}
64
66{
67}
68
70{
71 GlobalSysSolnType solvertype = m_locToGloMap.lock()->GetGlobalSysSolnType();
72 if (solvertype == eIterativeFull)
73 {
75 }
76 else if (solvertype == eIterativeStaticCond ||
77 solvertype == eIterativeMultiLevelStaticCond ||
78 solvertype == ePETScStaticCond ||
79 solvertype == ePETScMultiLevelStaticCond)
80 {
82 }
83 else
84 {
85 ASSERTL0(0, "Unsupported solver type");
86 }
87}
88
89/**
90 * Diagonal preconditioner computed by summing the relevant elements of
91 * the local matrix system.
92 */
94{
95 std::shared_ptr<MultiRegions::ExpList> expList =
96 ((m_linsys.lock())->GetLocMat()).lock();
97
98 auto asmMap = m_locToGloMap.lock();
99 int nGlobal = asmMap->GetNumGlobalCoeffs();
100 int nDir = asmMap->GetNumGlobalDirBndCoeffs();
101 int nInt = nGlobal - nDir;
102 int nElmt = expList->GetNumElmts();
103
104 Array<OneD, NekDouble> vOutput(nGlobal, 0.0);
105 for (int n = 0, cnt = 0; n < nElmt; ++n)
106 {
107 auto loc_mat = (m_linsys.lock())->GetBlock(n);
108 int loc_row = loc_mat->GetRows();
109 for (int i = 0; i < loc_row; ++i)
110 {
111 if (asmMap->GetLocalToGlobalSign(
112 cnt + i)) // check required for variable P
113 {
114 int gid1 = asmMap->GetLocalToGlobalMap(cnt + i);
115 vOutput[gid1] += (*loc_mat)(i, i);
116 }
117 }
118 cnt += loc_row;
119 }
120
121 // Assemble diagonal contributions across processes
122 asmMap->UniversalAssemble(vOutput);
123
125 Vmath::Sdiv(nInt, 1.0, &vOutput[nDir], 1, &m_diagonals[0], 1);
126}
127
128/**
129 * Diagonal preconditioner defined as the inverse of the main
130 * diagonal of the Schur complement
131 *
132 */
134{
135 auto asmMap = m_locToGloMap.lock();
136
137 int nGlobalBnd = asmMap->GetNumGlobalBndCoeffs();
138 int nDirBnd = asmMap->GetNumGlobalDirBndCoeffs();
139 int rows = nGlobalBnd - nDirBnd;
140
141 Array<OneD, NekDouble> vOutput(nGlobalBnd, 0.0);
142
143 // Extract diagonal contributions
145 for (unsigned int i = 0; i < rows; ++i)
146 {
147 vOutput[nDirBnd + i] = diagonals[i];
148 }
149
150 // Assemble diagonal contributions across processes
151 asmMap->UniversalAssembleBnd(vOutput);
152
154 Vmath::Sdiv(rows, 1.0, &vOutput[nDirBnd], 1, &m_diagonals[0], 1);
155}
156
157/**
158 *
159 */
161 const Array<OneD, NekDouble> &pInput, Array<OneD, NekDouble> &pOutput,
162 const bool &IsLocal)
163{
164 auto asmMap = m_locToGloMap.lock();
165
166 GlobalSysSolnType solvertype = asmMap->GetGlobalSysSolnType();
167
168 bool isFull = solvertype == eIterativeFull ? true : false;
169
170 int nGlobal = (isFull) ? asmMap->GetNumGlobalCoeffs()
171 : asmMap->GetNumGlobalBndCoeffs();
172 int nDir = asmMap->GetNumGlobalDirBndCoeffs();
173 int nNonDir = nGlobal - nDir;
174
175 if (IsLocal)
176 {
177 Array<OneD, NekDouble> wk(nGlobal);
178 (isFull) ? asmMap->Assemble(pInput, wk)
179 : asmMap->AssembleBnd(pInput, wk);
180 Vmath::Vmul(nNonDir, wk.data() + nDir, 1, m_diagonals.data(), 1,
181 wk.data() + nDir, 1);
182 Vmath::Zero(nDir, wk, 1);
183 (isFull) ? asmMap->GlobalToLocal(wk, pOutput)
184 : asmMap->GlobalToLocalBnd(wk, pOutput);
185 }
186 else
187 {
188 Vmath::Vmul(nNonDir, &pInput[0], 1, &m_diagonals[0], 1, &pOutput[0], 1);
189 }
190}
191
194 "Null", PreconditionerNull::create, "No Preconditioning");
195
196/**
197 * @class Null Preconditioner
198 *
199 * This class implements no preconditioning for the conjugate
200 * gradient matrix solver.
201 */
203 const std::shared_ptr<GlobalLinSys> &plinsys,
204 const AssemblyMapSharedPtr &pLocToGloMap)
205 : Preconditioner(plinsys, pLocToGloMap)
206{
207}
208
209/**
210 *
211 */
213{
214}
215
216/**
217 *
218 */
220{
221}
222
223/**
224 *
225 */
227 const Array<OneD, NekDouble> &pInput, Array<OneD, NekDouble> &pOutput,
228 const bool &isLocal)
229
230{
231 auto asmMap = m_locToGloMap.lock();
232
233 GlobalSysSolnType solvertype = asmMap->GetGlobalSysSolnType();
234
235 bool isFull = solvertype == eIterativeFull ? true : false;
236
237 int nGlobal = (isFull) ? asmMap->GetNumGlobalCoeffs()
238 : asmMap->GetNumGlobalBndCoeffs();
239 int nDir = asmMap->GetNumGlobalDirBndCoeffs();
240
241 if (isLocal)
242 {
243 Array<OneD, NekDouble> wk(nGlobal);
244 (isFull) ? asmMap->Assemble(pInput, wk)
245 : asmMap->AssembleBnd(pInput, wk);
246 Vmath::Zero(nDir, wk, 1);
247 (isFull) ? asmMap->GlobalToLocal(wk, pOutput)
248 : asmMap->GlobalToLocalBnd(wk, pOutput);
249 }
250 else
251 {
252 Vmath::Vcopy(pInput.size(), pInput, 1, pOutput, 1);
253 }
254}
255
256/**
257 * Registers the class with the Factory.
258 */
261 "Jacobi", PreconditionerJacobi::create, "Jacobi Preconditioning");
262/**
263 * @class PreconditionerJacobi
264 *
265 * This class implements jacobi preconditioning building on diagonal version
266 * above
267 */
268
270 const std::shared_ptr<GlobalLinSys> &plinsys,
271 const AssemblyMapSharedPtr &pLocToGloMap)
272 : PreconditionerDiagonal(plinsys, pLocToGloMap)
273{
274}
275
277{
278}
279
281{
283
284 auto expList = ((m_linsys.lock())->GetLocMat()).lock();
285 std::shared_ptr<LibUtilities::SessionReader> session =
286 expList->GetSession();
287
288 std::string var = m_locToGloMap.lock()->GetVariable();
289
290 if (session->DefinesGlobalSysSolnInfo(var, "JacobiIterations"))
291 {
292 m_niter = std::stoi(
293 session->GetGlobalSysSolnInfo(var, "JacobiIterations").c_str());
294 }
295 else
296 {
297 m_niter = 3;
298 }
299}
300
301/**
302 *
303 */
305 const Array<OneD, NekDouble> &pInput, Array<OneD, NekDouble> &pOutput,
306 const bool &IsLocal)
307{
308 auto asmMap = m_locToGloMap.lock();
309
310 GlobalSysSolnType solvertype = asmMap->GetGlobalSysSolnType();
311
312 int nGlobal = solvertype == eIterativeFull
313 ? asmMap->GetNumGlobalCoeffs()
314 : asmMap->GetNumGlobalBndCoeffs();
315 int nDir = asmMap->GetNumGlobalDirBndCoeffs();
316 int nNonDir = nGlobal - nDir;
317 Array<OneD, NekDouble> wk(nGlobal);
318
319 if (IsLocal)
320 {
321 int nLocal = solvertype == eIterativeFull
322 ? asmMap->GetNumLocalCoeffs()
323 : asmMap->GetNumLocalBndCoeffs();
324
325 Array<OneD, NekDouble> wk1(nLocal);
326 asmMap->Assemble(pInput, wk);
327 Vmath::Vmul(nNonDir, wk.data() + nDir, 1, m_diagonals.data(), 1,
328 wk.data() + nDir, 1);
329 Vmath::Zero(nDir, wk, 1);
330
331 for (int n = 1; n < m_niter; ++n)
332 {
333 asmMap->GlobalToLocal(wk, pOutput);
334
335 // do Ax operator
336 std::dynamic_pointer_cast<GlobalLinSysIterative>(m_linsys.lock())
337 ->DoMatrixMultiply(pOutput, wk1);
338
339 Vmath::Vsub(nLocal, pInput, 1, wk1, 1, wk1, 1);
340
341 asmMap->Assemble(wk1, pOutput);
342 Vmath::Vvtvp(nNonDir, pOutput.data() + nDir, 1, m_diagonals.data(),
343 1, wk.data() + nDir, 1, wk.data() + nDir, 1);
344 }
345
346 asmMap->GlobalToLocal(wk, pOutput);
347 }
348 else
349 {
350 Array<OneD, NekDouble> wk1(nGlobal);
351 Vmath::Vmul(nNonDir, pInput.data(), 1, m_diagonals.data(), 1,
352 wk.data() + nDir, 1);
353 Vmath::Zero(nDir, wk, 1);
354
355 for (int n = 1; n < m_niter; ++n)
356 {
357 // do Ax operator
358 std::dynamic_pointer_cast<GlobalLinSysIterative>(m_linsys.lock())
359 ->DoMatrixMultiply(wk, wk1);
360
361 // b - Ax
362 Vmath::Vsub(nNonDir, pInput.data(), 1, wk1.data() + nDir, 1,
363 wk1.data() + nDir, 1);
364
365 // new sol = 1/diag (b-Ax) + old sol
366 Vmath::Vvtvp(nNonDir, wk1.data() + nDir, 1, m_diagonals.data(), 1,
367 wk.data() + nDir, 1, wk.data() + nDir, 1);
368 }
369
370 Vmath::Vcopy(nNonDir, wk.data() + nDir, 1, pOutput.data(), 1);
371 }
372}
373
374} // namespace Nektar::MultiRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
PreconditionerDiagonal(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &IsLocal=false) override
Apply a preconditioner to the conjugate gradient method.
const std::weak_ptr< GlobalLinSys > m_linsys
Array< OneD, NekDouble > AssembleStaticCondGlobalDiagonals()
Performs global assembly of diagonal entries to global Schur complement matrix.
std::weak_ptr< AssemblyMap > m_locToGloMap
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
PreconditionerJacobi(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
static std::string className
Name of class.
void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &IsLocal=false) override
Apply a preconditioner to the conjugate gradient method.
void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false) override
Apply a preconditioner to the conjugate gradient method.
static std::string className
Name of class.
PreconditionerNull(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
PreconFactory & GetPreconFactory()
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:50
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
Definition: Vmath.hpp:154
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
STL namespace.