Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  solvertype == ePETScStaticCond ||
84  solvertype == ePETScMultiLevelStaticCond)
85  {
87  }
88  else
89  {
90  ASSERTL0(0,"Unsupported solver type");
91  }
92  }
93 
94  /**
95  * Diagonal preconditioner computed by summing the relevant elements of
96  * the local matrix system.
97  */
99  {
100  boost::shared_ptr<MultiRegions::ExpList> expList =
101  ((m_linsys.lock())->GetLocMat()).lock();
102 
104 
105  int i,j,n,cnt,gid1,gid2;
106  NekDouble sign1,sign2,value;
107  int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
108  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
109  int nInt = nGlobal - nDir;
110 
111  // fill global matrix
112  DNekScalMatSharedPtr loc_mat;
113  Array<OneD, NekDouble> vOutput(nGlobal,0.0);
114 
115  int loc_row;
116  int nElmt = expList->GetNumElmts();
117  for(n = cnt = 0; n < nElmt; ++n)
118  {
119  loc_mat = (m_linsys.lock())->GetBlock(expList->GetOffset_Elmt_Id(n));
120  loc_row = loc_mat->GetRows();
121 
122  for(i = 0; i < loc_row; ++i)
123  {
124  gid1 = m_locToGloMap->GetLocalToGlobalMap(cnt + i) - nDir;
125  sign1 = m_locToGloMap->GetLocalToGlobalSign(cnt + i);
126  if(gid1 >= 0)
127  {
128  for(j = 0; j < loc_row; ++j)
129  {
130  gid2 = m_locToGloMap->GetLocalToGlobalMap(cnt + j)
131  - nDir;
132  sign2 = m_locToGloMap->GetLocalToGlobalSign(cnt + j);
133  if(gid2 == gid1)
134  {
135  // When global matrix is symmetric,
136  // only add the value for the upper
137  // triangular part in order to avoid
138  // entries to be entered twice
139  value = vOutput[gid1 + nDir]
140  + sign1*sign2*(*loc_mat)(i,j);
141  vOutput[gid1 + nDir] = value;
142  }
143  }
144  }
145  }
146  cnt += loc_row;
147  }
148 
149  // Assemble diagonal contributions across processes
150  m_locToGloMap->UniversalAssemble(vOutput);
151 
153  Vmath::Sdiv(nInt, 1.0, &vOutput[nDir], 1, &m_diagonals[0], 1);
154  }
155 
156  /**
157  * Diagonal preconditioner defined as the inverse of the main
158  * diagonal of the Schur complement
159  *
160  */
162  {
163  int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
164  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
165  int rows = nGlobalBnd - nDirBnd;
166 
167  Array<OneD, NekDouble> vOutput(nGlobalBnd,0.0);
168 
169  // Extract diagonal contributions
171  for (unsigned int i = 0; i < rows; ++i)
172  {
173  vOutput[nDirBnd + i] = diagonals[i];
174  }
175 
176  // Assemble diagonal contributions across processes
177  m_locToGloMap->UniversalAssembleBnd(vOutput);
178 
180  Vmath::Sdiv(rows, 1.0, &vOutput[nDirBnd], 1, &m_diagonals[0], 1);
181  }
182 
183  /**
184  *
185  */
187  const Array<OneD, NekDouble>& pInput,
188  Array<OneD, NekDouble>& pOutput)
189  {
190  GlobalSysSolnType solvertype =
191  m_locToGloMap->GetGlobalSysSolnType();
192 
193  int nGlobal = solvertype == eIterativeFull ?
194  m_locToGloMap->GetNumGlobalCoeffs() :
195  m_locToGloMap->GetNumGlobalBndCoeffs();
196  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
197  int nNonDir = nGlobal-nDir;
198  Vmath::Vmul(nNonDir, &pInput[0], 1, &m_diagonals[0], 1, &pOutput[0], 1);
199  }
200 
203  "Null",
205  "No Preconditioning");
206 
207  /**
208  * @class Null Preconditioner
209  *
210  * This class implements no preconditioning for the conjugate
211  * gradient matrix solver.
212  */
214  const boost::shared_ptr<GlobalLinSys> &plinsys,
215  const AssemblyMapSharedPtr &pLocToGloMap)
216  : Preconditioner(plinsys, pLocToGloMap),
217  m_preconType(pLocToGloMap->GetPreconType())
218  {
219  }
220 
221  /**
222  *
223  */
225  {
226  }
227 
228  /**
229  *
230  */
232  {
233  }
234 
235  /**
236  *
237  */
239  const Array<OneD, NekDouble>& pInput,
240  Array<OneD, NekDouble>& pOutput)
241  {
242  Vmath::Vcopy(pInput.num_elements(), pInput, 1, pOutput, 1);
243  }
244 
245 
246  }
247 }
248 
249 
250 
251 
252 
253 
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:1047
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