Nektar++
PreconditionerDiagonal.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Preconditioner.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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
40 #include <math.h>
41 
42 namespace Nektar
43 {
44  namespace MultiRegions
45  {
46  /**
47  * Registers the class with the Factory.
48  */
51  "Diagonal",
53  "Diagonal Preconditioning");
54  /**
55  * @class Preconditioner
56  *
57  * This class implements diagonal preconditioning for the conjugate
58  * gradient matrix solver.
59  */
60 
62  const boost::shared_ptr<GlobalLinSys> &plinsys,
63  const AssemblyMapSharedPtr &pLocToGloMap)
64  : Preconditioner(plinsys, pLocToGloMap),
65  m_preconType(pLocToGloMap->GetPreconType())
66  {
67  }
68 
70  {
71  }
72 
74  {
75  GlobalSysSolnType solvertype =
76  m_locToGloMap->GetGlobalSysSolnType();
77  if (solvertype == eIterativeFull)
78  {
80  }
81  else if(solvertype == eIterativeStaticCond ||
82  solvertype == eIterativeMultiLevelStaticCond)
83  {
85  }
86  else
87  {
88  ASSERTL0(0,"Unsupported solver type");
89  }
90  }
91 
92  /**
93  * Diagonal preconditioner computed by summing the relevant elements of
94  * the local matrix system.
95  */
97  {
98  boost::shared_ptr<MultiRegions::ExpList> expList =
99  ((m_linsys.lock())->GetLocMat()).lock();
100 
102 
103  int i,j,n,cnt,gid1,gid2;
104  NekDouble sign1,sign2,value;
105  int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
106  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
107  int nInt = nGlobal - nDir;
108 
109  // fill global matrix
110  DNekScalMatSharedPtr loc_mat;
111  Array<OneD, NekDouble> vOutput(nGlobal,0.0);
112 
113  int loc_row;
114  int nElmt = expList->GetNumElmts();
115  for(n = cnt = 0; n < nElmt; ++n)
116  {
117  loc_mat = (m_linsys.lock())->GetBlock(expList->GetOffset_Elmt_Id(n));
118  loc_row = loc_mat->GetRows();
119 
120  for(i = 0; i < loc_row; ++i)
121  {
122  gid1 = m_locToGloMap->GetLocalToGlobalMap(cnt + i) - nDir;
123  sign1 = m_locToGloMap->GetLocalToGlobalSign(cnt + i);
124  if(gid1 >= 0)
125  {
126  for(j = 0; j < loc_row; ++j)
127  {
128  gid2 = m_locToGloMap->GetLocalToGlobalMap(cnt + j)
129  - nDir;
130  sign2 = m_locToGloMap->GetLocalToGlobalSign(cnt + j);
131  if(gid2 == gid1)
132  {
133  // When global matrix is symmetric,
134  // only add the value for the upper
135  // triangular part in order to avoid
136  // entries to be entered twice
137  value = vOutput[gid1 + nDir]
138  + sign1*sign2*(*loc_mat)(i,j);
139  vOutput[gid1 + nDir] = value;
140  }
141  }
142  }
143  }
144  cnt += loc_row;
145  }
146 
147  // Assemble diagonal contributions across processes
148  m_locToGloMap->UniversalAssemble(vOutput);
149 
151  Vmath::Sdiv(nInt, 1.0, &vOutput[nDir], 1, &m_diagonals[0], 1);
152  }
153 
154  /**
155  * Diagonal preconditioner defined as the inverse of the main
156  * diagonal of the Schur complement
157  *
158  */
160  {
161  int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
162  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
163  int rows = nGlobalBnd - nDirBnd;
164 
165  Array<OneD, NekDouble> vOutput(nGlobalBnd,0.0);
166 
167  // Extract diagonal contributions
169  for (unsigned int i = 0; i < rows; ++i)
170  {
171  vOutput[nDirBnd + i] = diagonals[i];
172  }
173 
174  // Assemble diagonal contributions across processes
175  m_locToGloMap->UniversalAssembleBnd(vOutput);
176 
178  Vmath::Sdiv(rows, 1.0, &vOutput[nDirBnd], 1, &m_diagonals[0], 1);
179  }
180 
181  /**
182  *
183  */
185  const Array<OneD, NekDouble>& pInput,
186  Array<OneD, NekDouble>& pOutput)
187  {
188  GlobalSysSolnType solvertype =
189  m_locToGloMap->GetGlobalSysSolnType();
190 
191  int nGlobal = solvertype == eIterativeFull ?
192  m_locToGloMap->GetNumGlobalCoeffs() :
193  m_locToGloMap->GetNumGlobalBndCoeffs();
194  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
195  int nNonDir = nGlobal-nDir;
196  Vmath::Vmul(nNonDir, &pInput[0], 1, &m_diagonals[0], 1, &pOutput[0], 1);
197  }
198 
201  "Null",
203  "No Preconditioning");
204 
205  /**
206  * @class Null Preconditioner
207  *
208  * This class implements no preconditioning for the conjugate
209  * gradient matrix solver.
210  */
212  const boost::shared_ptr<GlobalLinSys> &plinsys,
213  const AssemblyMapSharedPtr &pLocToGloMap)
214  : Preconditioner(plinsys, pLocToGloMap),
215  m_preconType(pLocToGloMap->GetPreconType())
216  {
217  }
218 
219  /**
220  *
221  */
223  {
224  }
225 
226  /**
227  *
228  */
230  {
231  }
232 
233  /**
234  *
235  */
237  const Array<OneD, NekDouble>& pInput,
238  Array<OneD, NekDouble>& pOutput)
239  {
240  Vmath::Vcopy(pInput.num_elements(), pInput, 1, pOutput, 1);
241  }
242 
243 
244  }
245 }
246 
247 
248 
249 
250 
251 
boost::shared_ptr< AssemblyMap > m_locToGloMap
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Apply a preconditioner to the conjugate gradient method.
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Apply a preconditioner to the conjugate gradient method.
PreconFactory & GetPreconFactory()
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
const boost::weak_ptr< GlobalLinSys > m_linsys
Array< OneD, NekDouble > AssembleStaticCondGlobalDiagonals()
Performs global assembly of diagonal entries to global Schur complement matrix.
double NekDouble
PreconditionerDiagonal(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
PreconditionerNull(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
static PreconditionerSharedPtr create(const boost::shared_ptr< GlobalLinSys > &plinsys, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
static std::string className
Name of class.
static PreconditionerSharedPtr create(const boost::shared_ptr< GlobalLinSys > &plinsys, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
static std::string className
Name of class.
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.cpp:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215