Nektar++
GlobalLinSysIterative.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: GlobalLinSysIterative.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: GlobalLinSysIterative definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40  namespace MultiRegions
41  {
42  /**
43  * @class GlobalLinSysIterative
44  *
45  * Solves a linear system using iterative methods.
46  */
47 
48  /// Constructor for full direct matrix solve.
50  const GlobalLinSysKey &pKey,
51  const boost::weak_ptr<ExpList> &pExpList,
52  const boost::shared_ptr<AssemblyMap>
53  &pLocToGloMap)
54  : GlobalLinSys(pKey, pExpList, pLocToGloMap),
55  m_rhs_magnitude(NekConstants::kNekUnsetDouble),
57  m_totalIterations(0),
58  m_useProjection(false),
59  m_numPrevSols(0)
60  {
62  = pExpList.lock()->GetSession();
63 
64  m_tolerance = pLocToGloMap->GetIterativeTolerance();
65 
66  LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
67  m_root = (vComm->GetRank())? false : true;
68  m_verbose = (vSession->DefinesCmdLineArgument("verbose"))? true :false;
69 
70  int successiveRHS;
71 
72  if((successiveRHS = pLocToGloMap->GetSuccessiveRHS()))
73  {
74  m_prevLinSol.set_capacity(successiveRHS);
75  m_useProjection = true;
76  }
77  else
78  {
79  m_useProjection = false;
80  }
81  }
82 
84  {
85  }
86 
87 
88  /**
89  *
90  */
92  const int nGlobal,
93  const Array<OneD,const NekDouble> &pInput,
94  Array<OneD, NekDouble> &pOutput,
95  const AssemblyMapSharedPtr &plocToGloMap,
96  const int nDir)
97  {
98  if (m_useProjection)
99  {
100  DoAconjugateProjection(nGlobal, pInput, pOutput, plocToGloMap, nDir);
101  }
102  else
103  {
104  // applying plain Conjugate Gradient
105  DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
106  }
107  }
108 
109 
110  /**
111  * This method implements A-conjugate projection technique
112  * in order to speed up successive linear solves with
113  * right-hand sides arising from time-dependent discretisations.
114  * (P.F.Fischer, Comput. Methods Appl. Mech. Engrg. 163, 1998)
115  */
117  const int nGlobal,
118  const Array<OneD,const NekDouble> &pInput,
119  Array<OneD, NekDouble> &pOutput,
120  const AssemblyMapSharedPtr &plocToGloMap,
121  const int nDir)
122  {
123  // Get the communicator for performing data exchanges
125  = m_expList.lock()->GetComm()->GetRowComm();
126 
127  // Get vector sizes
128  int nNonDir = nGlobal - nDir;
130 
131  if (0 == m_numPrevSols)
132  {
133  // no previous solutions found, call CG
134 
135  DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
136 
137  UpdateKnownSolutions(nGlobal, pOutput, nDir);
138  }
139  else
140  {
141  // Create NekVector wrappers for linear algebra operations
142  NekVector<NekDouble> b (nNonDir, pInput + nDir, eWrapper);
143  NekVector<NekDouble> x (nNonDir, tmp = pOutput + nDir, eWrapper);
144 
145  // check the input vector (rhs) is not zero
146 
147  NekDouble rhsNorm = Vmath::Dot2(nNonDir,
148  pInput + nDir,
149  pInput + nDir,
150  m_map + nDir);
151 
152  vComm->AllReduce(rhsNorm, Nektar::LibUtilities::ReduceSum);
153 
154  if (rhsNorm < NekConstants::kNekZeroTol)
155  {
156  Array<OneD, NekDouble> tmp = pOutput+nDir;
157  Vmath::Zero(nNonDir, tmp, 1);
158  return;
159  }
160 
161  // Allocate array storage
162  Array<OneD, NekDouble> px_s (nGlobal, 0.0);
163  Array<OneD, NekDouble> pb_s (nGlobal, 0.0);
164  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
165  Array<OneD, NekDouble> tmpx_s (nGlobal, 0.0);
166 
167  NekVector<NekDouble> pb (nNonDir, tmp = pb_s + nDir, eWrapper);
168  NekVector<NekDouble> px (nNonDir, tmp = px_s + nDir, eWrapper);
169  NekVector<NekDouble> tmpAx (nNonDir, tmp = tmpAx_s + nDir, eWrapper);
170  NekVector<NekDouble> tmpx (nNonDir, tmp = tmpx_s + nDir, eWrapper);
171 
172 
173  // notation follows the paper cited:
174  // \alpha_i = \tilda{x_i}^T b^n
175  // projected x, px = \sum \alpha_i \tilda{x_i}
176 
177  Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
178  for (int i = 0; i < m_prevLinSol.size(); i++)
179  {
180  alpha[i] = Vmath::Dot2(nNonDir,
181  m_prevLinSol[i],
182  pInput + nDir,
183  m_map + nDir);
184  }
185  vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
186 
187  for (int i = 0; i < m_prevLinSol.size(); i++)
188  {
189  if (alpha[i] < NekConstants::kNekZeroTol)
190  {
191  continue;
192  }
193 
194  NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
195  px += alpha[i] * xi;
196  }
197 
198  // pb = b^n - A px
199  Vmath::Vcopy(nNonDir,
200  pInput.get() + nDir, 1,
201  pb_s.get() + nDir, 1);
202 
203  v_DoMatrixMultiply(px_s, tmpAx_s);
204 
205  pb -= tmpAx;
206 
207 
208  // solve the system with projected rhs
209  DoConjugateGradient(nGlobal, pb_s, tmpx_s, plocToGloMap, nDir);
210 
211 
212  // remainder solution + projection of previous solutions
213  x = tmpx + px;
214 
215  // save the auxiliary solution to prev. known solutions
216  UpdateKnownSolutions(nGlobal, tmpx_s, nDir);
217  }
218  }
219 
220 
221  /**
222  * Calculating A-norm of an input vector,
223  * A-norm(x) := sqrt( < x, Ax > )
224  */
226  const int nGlobal,
227  const Array<OneD,const NekDouble> &in,
228  const int nDir)
229  {
230  // Get the communicator for performing data exchanges
232  = m_expList.lock()->GetComm()->GetRowComm();
233 
234  // Get vector sizes
235  int nNonDir = nGlobal - nDir;
236 
237  // Allocate array storage
238  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
239 
240  v_DoMatrixMultiply(in, tmpAx_s);
241 
242  NekDouble anorm_sq = Vmath::Dot2(nNonDir,
243  in + nDir,
244  tmpAx_s + nDir,
245  m_map + nDir);
246  vComm->AllReduce(anorm_sq, Nektar::LibUtilities::ReduceSum);
247  return std::sqrt(anorm_sq);
248  }
249 
250  /**
251  * Updates the storage of previously known solutions.
252  * Performs normalisation of input vector wrt A-norm.
253  */
255  const int nGlobal,
256  const Array<OneD,const NekDouble> &newX,
257  const int nDir)
258  {
259  // Get vector sizes
260  int nNonDir = nGlobal - nDir;
261 
262  // Get the communicator for performing data exchanges
264  = m_expList.lock()->GetComm()->GetRowComm();
265 
266  // Check the solution is non-zero
267  NekDouble solNorm = Vmath::Dot2(nNonDir,
268  newX + nDir,
269  newX + nDir,
270  m_map + nDir);
271  vComm->AllReduce(solNorm, Nektar::LibUtilities::ReduceSum);
272 
273  if (solNorm < NekConstants::kNekZeroTol)
274  {
275  return;
276  }
277 
278 
279  // Allocate array storage
280  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
281  Array<OneD, NekDouble> px_s (nGlobal, 0.0);
282  Array<OneD, NekDouble> tmp1, tmp2;
283 
284  // Create NekVector wrappers for linear algebra operations
285  NekVector<NekDouble> px (nNonDir, tmp1 = px_s + nDir, eWrapper);
286  NekVector<NekDouble> tmpAx (nNonDir, tmp2 = tmpAx_s + nDir, eWrapper);
287 
288 
289  // calculating \tilda{x} - sum \alpha_i\tilda{x}_i
290 
291  Vmath::Vcopy(nNonDir,
292  tmp1 = newX + nDir, 1,
293  tmp2 = px_s + nDir, 1);
294 
295  if (m_prevLinSol.size() > 0)
296  {
297  v_DoMatrixMultiply(newX, tmpAx_s);
298  }
299 
300  Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
301  for (int i = 0; i < m_prevLinSol.size(); i++)
302  {
303  alpha[i] = Vmath::Dot2(nNonDir,
304  m_prevLinSol[i],
305  tmpAx_s + nDir,
306  m_map + nDir);
307  }
308  vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
309 
310  for (int i = 0; i < m_prevLinSol.size(); i++)
311  {
312  if (alpha[i] < NekConstants::kNekZeroTol)
313  {
314  continue;
315  }
316 
317  NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
318  px -= alpha[i] * xi;
319  }
320 
321 
322  // Some solutions generated by CG are identical zeros, see
323  // solutions generated for Test_Tet_equitri.xml (IncNavierStokesSolver).
324  // Not going to store identically zero solutions.
325 
326  NekDouble anorm = CalculateAnorm(nGlobal, px_s, nDir);
327  if (anorm < NekConstants::kNekZeroTol)
328  {
329  return;
330  }
331 
332  // normalisation of new solution
333  Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
334 
335  // updating storage with non-Dirichlet-dof part of new solution vector
336  m_prevLinSol.push_back(px_s + nDir);
337  m_numPrevSols++;
338  }
339 
340 
341 
342  /**  
343  * Solve a global linear system using the conjugate gradient method.  
344  * We solve only for the non-Dirichlet modes. The operator is evaluated  
345  * using an auxiliary function v_DoMatrixMultiply defined by the  
346  * specific solver. Distributed math routines are used to support  
347  * parallel execution of the solver.  
348  *  
349  * The implemented algorithm uses a reduced-communication reordering of  
350  * the standard PCG method (Demmel, Heath and Vorst, 1993)  
351  *  
352  * @param pInput Input residual of all DOFs.  
353  * @param pOutput Solution vector of all DOFs.  
354  */
356  const int nGlobal,
357  const Array<OneD,const NekDouble> &pInput,
358  Array<OneD, NekDouble> &pOutput,
359  const AssemblyMapSharedPtr &plocToGloMap,
360  const int nDir)
361  {
362  if (!m_precon)
363  {
365  = plocToGloMap->GetPreconType();
366  std::string PreconType
368  v_UniqueMap();
370  PreconType,GetSharedThisPtr(),plocToGloMap);
371  m_precon->BuildPreconditioner();
372  }
373 
374  // Get the communicator for performing data exchanges
376  = m_expList.lock()->GetComm()->GetRowComm();
377 
378  // Get vector sizes
379  int nNonDir = nGlobal - nDir;
380 
381  // Allocate array storage
382  Array<OneD, NekDouble> w_A (nGlobal, 0.0);
383  Array<OneD, NekDouble> s_A (nGlobal, 0.0);
384  Array<OneD, NekDouble> p_A (nNonDir, 0.0);
385  Array<OneD, NekDouble> r_A (nNonDir, 0.0);
386  Array<OneD, NekDouble> q_A (nNonDir, 0.0);
388 
389  // Create NekVector wrappers for linear algebra operations
390  NekVector<NekDouble> in (nNonDir,pInput + nDir, eWrapper);
391  NekVector<NekDouble> out(nNonDir,tmp = pOutput + nDir,eWrapper);
392  NekVector<NekDouble> w (nNonDir,tmp = w_A + nDir, eWrapper);
393  NekVector<NekDouble> s (nNonDir,tmp = s_A + nDir, eWrapper);
394  NekVector<NekDouble> p (nNonDir,p_A, eWrapper);
395  NekVector<NekDouble> r (nNonDir,r_A, eWrapper);
396  NekVector<NekDouble> q (nNonDir,q_A, eWrapper);
397 
398  int k;
399  NekDouble alpha, beta, rho, rho_new, mu, eps, min_resid;
400  Array<OneD, NekDouble> vExchange(3,0.0);
401 
402  // Copy initial residual from input
403  r = in;
404  // zero homogeneous out array ready for solution updates
405  // Should not be earlier in case input vector is same as
406  // output and above copy has been peformed
407  Vmath::Zero(nNonDir,tmp = pOutput + nDir,1);
408 
409 
410  // evaluate initial residual error for exit check
411  vExchange[2] = Vmath::Dot2(nNonDir,
412  r_A,
413  r_A,
414  m_map + nDir);
415 
416  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
417 
418  eps = vExchange[2];
419 
421  {
422  m_rhs_magnitude = 1.0/vExchange[2];
423  }
424 
425  // If input residual is less than tolerance skip solve.
427  {
428  if (m_verbose && m_root)
429  {
430  cout << "CG iterations made = " << m_totalIterations
431  << " using tolerance of " << m_tolerance
432  << " (error = " << sqrt(eps/m_rhs_magnitude) << ")" << endl;
433  }
434  m_rhs_magnitude = NekConstants::kNekUnsetDouble;
435  return;
436  }
437 
438  m_totalIterations = 1;
439  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
440 
441  v_DoMatrixMultiply(w_A, s_A);
442 
443  k = 0;
444 
445  vExchange[0] = Vmath::Dot2(nNonDir,
446  r_A,
447  w_A + nDir,
448  m_map + nDir);
449 
450  vExchange[1] = Vmath::Dot2(nNonDir,
451  s_A + nDir,
452  w_A + nDir,
453  m_map + nDir);
454 
455  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
456 
457  rho = vExchange[0];
458  mu = vExchange[1];
459  min_resid = m_rhs_magnitude;
460  beta = 0.0;
461  alpha = rho/mu;
462 
463  // Continue until convergence
464  while (true)
465  {
466  ASSERTL0(k < 5000,
467  "Exceeded maximum number of iterations (5000)");
468 
469  // Compute new search direction p_k, q_k
470  Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
471  Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
472 
473  // Update solution x_{k+1}
474  Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
475 
476  // Update residual vector r_{k+1}
477  Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
478 
479  // Apply preconditioner
480  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
481 
482  // Perform the method-specific matrix-vector multiply operation.
483  v_DoMatrixMultiply(w_A, s_A);
484 
485  // <r_{k+1}, w_{k+1}>
486  vExchange[0] = Vmath::Dot2(nNonDir,
487  r_A,
488  w_A + nDir,
489  m_map + nDir);
490  // <s_{k+1}, w_{k+1}>
491  vExchange[1] = Vmath::Dot2(nNonDir,
492  s_A + nDir,
493  w_A + nDir,
494  m_map + nDir);
495 
496  // <r_{k+1}, r_{k+1}>
497  vExchange[2] = Vmath::Dot2(nNonDir,
498  r_A,
499  r_A,
500  m_map + nDir);
501 
502  // Perform inner-product exchanges
503  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
504 
505  rho_new = vExchange[0];
506  mu = vExchange[1];
507  eps = vExchange[2];
508 
510  // test if norm is within tolerance
511  if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
512  {
513  if (m_verbose && m_root)
514  {
515  cout << "CG iterations made = " << m_totalIterations
516  << " using tolerance of " << m_tolerance
517  << " (error = " << sqrt(eps/m_rhs_magnitude) << ")"
518  //<< " (bb_inv = " << sqrt(1.0/m_rhs_magnitude) << ")"
519  << endl;
520  }
521  m_rhs_magnitude = NekConstants::kNekUnsetDouble;
522  break;
523  }
524  min_resid = min(min_resid, eps);
525 
526  // Compute search direction and solution coefficients
527  beta = rho_new/rho;
528  alpha = rho_new/(mu - rho_new*beta/alpha);
529  rho = rho_new;
530  k++;
531  }
532  }
533 
535  {
536 
537  Array<OneD, NekDouble> vExchange(1);
538  vExchange[0] = Vmath::Dot(pIn.GetDimension(),&pIn[0],&pIn[0]);
539 
540  m_expList.lock()->GetComm()->GetRowComm()->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
541  m_rhs_magnitude = (vExchange[0] > 1e-6)? vExchange[0]:1.0;
542  }
543 
544  }
545 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
void Set_Rhs_Magnitude(const NekVector< NekDouble > &pIn)
void UpdateKnownSolutions(const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs)
bool m_root
Provide verbose output and root if parallel.
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
bool m_useProjection
Whether to apply projection technique.
PreconFactory & GetPreconFactory()
boost::circular_buffer< Array< OneD, NekDouble > > m_prevLinSol
Storage for solutions to previous linear problems.
virtual void v_DoMatrixMultiply(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
NekDouble CalculateAnorm(const int nGlobal, const Array< OneD, const NekDouble > &in, const int nDir)
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
static const NekDouble kNekZeroTol
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
static PreconditionerSharedPtr NullPreconditionerSharedPtr
NekDouble m_tolerance
Tolerance of iterative solver.
int m_numPrevSols
Total counter of previous solutions.
Array< OneD, int > m_map
Global to universal unique map.
boost::shared_ptr< GlobalLinSys > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: GlobalLinSys.h:103
double NekDouble
GlobalLinSysIterative(const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
static const NekDouble kNekUnsetDouble
Describe a linear system.
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:891
unsigned int GetDimension() const
Returns the number of dimensions for the point.
Definition: NekVector.cpp:212
void DoAconjugateProjection(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
A-conjugate projection technique.
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
A global linear system.
Definition: GlobalLinSys.h:70
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:931
virtual void v_SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
Solve the matrix system.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
const char *const PreconditionerTypeMap[]
void DoConjugateGradient(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
Actual iterative solve.
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:125