Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
48  /**
49  * Registers the class with the Factory.
50  */
51  string PreconditionerDiagonal::className
53  "Diagonal",
54  PreconditionerDiagonal::create,
55  "Diagonal Preconditioning");
56  /**
57  * @class Preconditioner
58  *
59  * This class implements diagonal preconditioning for the conjugate
60  * gradient matrix solver.
61  */
62 
63  PreconditionerDiagonal::PreconditionerDiagonal(
64  const boost::shared_ptr<GlobalLinSys> &plinsys,
65  const AssemblyMapSharedPtr &pLocToGloMap)
66  : Preconditioner(plinsys, pLocToGloMap),
67  m_preconType(pLocToGloMap->GetPreconType())
68  {
69  }
70 
72  {
73  }
74 
76  {
77  GlobalSysSolnType solvertype =
78  m_locToGloMap->GetGlobalSysSolnType();
79  if (solvertype == eIterativeFull)
80  {
82  }
83  else if(solvertype == eIterativeStaticCond ||
84  solvertype == eIterativeMultiLevelStaticCond ||
85  solvertype == ePETScStaticCond ||
86  solvertype == ePETScMultiLevelStaticCond)
87  {
89  }
90  else
91  {
92  ASSERTL0(0,"Unsupported solver type");
93  }
94  }
95 
96  /**
97  * Diagonal preconditioner computed by summing the relevant elements of
98  * the local matrix system.
99  */
101  {
102  boost::shared_ptr<MultiRegions::ExpList> expList =
103  ((m_linsys.lock())->GetLocMat()).lock();
104 
106 
107  int i,j,n,cnt,gid1,gid2;
108  NekDouble sign1,sign2,value;
109  int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
110  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
111  int nInt = nGlobal - nDir;
112 
113  // fill global matrix
114  DNekScalMatSharedPtr loc_mat;
115  Array<OneD, NekDouble> vOutput(nGlobal,0.0);
116 
117  int loc_row;
118  int nElmt = expList->GetNumElmts();
119  for(n = cnt = 0; n < nElmt; ++n)
120  {
121  loc_mat = (m_linsys.lock())->GetBlock(n);
122  loc_row = loc_mat->GetRows();
123 
124  for(i = 0; i < loc_row; ++i)
125  {
126  gid1 = m_locToGloMap->GetLocalToGlobalMap(cnt + i) - nDir;
127  sign1 = m_locToGloMap->GetLocalToGlobalSign(cnt + i);
128  if(gid1 >= 0)
129  {
130  for(j = 0; j < loc_row; ++j)
131  {
132  gid2 = m_locToGloMap->GetLocalToGlobalMap(cnt + j)
133  - nDir;
134  sign2 = m_locToGloMap->GetLocalToGlobalSign(cnt + j);
135  if(gid2 == gid1)
136  {
137  // When global matrix is symmetric,
138  // only add the value for the upper
139  // triangular part in order to avoid
140  // entries to be entered twice
141  value = vOutput[gid1 + nDir]
142  + sign1*sign2*(*loc_mat)(i,j);
143  vOutput[gid1 + nDir] = value;
144  }
145  }
146  }
147  }
148  cnt += loc_row;
149  }
150 
151  // Assemble diagonal contributions across processes
152  m_locToGloMap->UniversalAssemble(vOutput);
153 
155  Vmath::Sdiv(nInt, 1.0, &vOutput[nDir], 1, &m_diagonals[0], 1);
156  }
157 
158  /**
159  * Diagonal preconditioner defined as the inverse of the main
160  * diagonal of the Schur complement
161  *
162  */
164  {
165  int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
166  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
167  int rows = nGlobalBnd - nDirBnd;
168 
169  Array<OneD, NekDouble> vOutput(nGlobalBnd,0.0);
170 
171  // Extract diagonal contributions
173  for (unsigned int i = 0; i < rows; ++i)
174  {
175  vOutput[nDirBnd + i] = diagonals[i];
176  }
177 
178  // Assemble diagonal contributions across processes
179  m_locToGloMap->UniversalAssembleBnd(vOutput);
180 
182  Vmath::Sdiv(rows, 1.0, &vOutput[nDirBnd], 1, &m_diagonals[0], 1);
183  }
184 
185  /**
186  *
187  */
189  const Array<OneD, NekDouble>& pInput,
190  Array<OneD, NekDouble>& pOutput)
191  {
192  GlobalSysSolnType solvertype =
193  m_locToGloMap->GetGlobalSysSolnType();
194 
195  int nGlobal = solvertype == eIterativeFull ?
196  m_locToGloMap->GetNumGlobalCoeffs() :
197  m_locToGloMap->GetNumGlobalBndCoeffs();
198  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
199  int nNonDir = nGlobal-nDir;
200  Vmath::Vmul(nNonDir, &pInput[0], 1, &m_diagonals[0], 1, &pOutput[0], 1);
201  }
202 
205  "Null",
207  "No Preconditioning");
208 
209  /**
210  * @class Null Preconditioner
211  *
212  * This class implements no preconditioning for the conjugate
213  * gradient matrix solver.
214  */
216  const boost::shared_ptr<GlobalLinSys> &plinsys,
217  const AssemblyMapSharedPtr &pLocToGloMap)
218  : Preconditioner(plinsys, pLocToGloMap),
219  m_preconType(pLocToGloMap->GetPreconType())
220  {
221  }
222 
223  /**
224  *
225  */
227  {
228  }
229 
230  /**
231  *
232  */
234  {
235  }
236 
237  /**
238  *
239  */
241  const Array<OneD, NekDouble>& pInput,
242  Array<OneD, NekDouble>& pOutput)
243  {
244  Vmath::Vcopy(pInput.num_elements(), pInput, 1, pOutput, 1);
245  }
246 
247 
248  }
249 }
250 
251 
252 
253 
254 
255 
boost::shared_ptr< AssemblyMap > m_locToGloMap
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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()
STL namespace.
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:271
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
PreconditionerNull(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
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:1061
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:183
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215