Nektar++
PreconditionerBlock.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PreconditionerBlock.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: Block Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
42 #include <LocalRegions/MatrixKey.h>
43 #include <LocalRegions/SegExp.h>
44 #include <math.h>
45 
46 namespace Nektar
47 {
48  namespace MultiRegions
49  {
50  /**
51  * Registers the class with the Factory.
52  */
55  "Block",
57  "Block Diagonal Preconditioning");
58  /**
59  * @class Block Preconditioner
60  *
61  * This class implements Block preconditioning for the conjugate
62  * gradient matrix solver.
63  */
64 
66  const boost::shared_ptr<GlobalLinSys> &plinsys,
67  const AssemblyMapSharedPtr &pLocToGloMap)
68  : Preconditioner(plinsys, pLocToGloMap),
69  m_linsys(plinsys),
70  m_preconType(pLocToGloMap->GetPreconType()),
71  m_locToGloMap(pLocToGloMap)
72  {
73  }
74 
76  {
77  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
79  "Solver type not valid");
80  }
81 
83  {
84  GlobalLinSysKey key = m_linsys.lock()->GetKey();
85 
86  // Different setup for HDG and CG.
88  {
90  }
91  else
92  {
94  }
95  }
96 
97  /**
98  * \brief Construct a block preconditioner from \f$\mathbf{S}_{1}\f$ for
99  * the continuous Galerkin system.
100  *
101  * The preconditioner is defined as:
102  *
103  * \f[\mathbf{M}^{-1}=\left[\begin{array}{ccc}
104  * Diag[(\mathbf{S_{1}})_{vv}] & & \\ & (\mathbf{S}_{1})_{eb} & \\ & &
105  * (\mathbf{S}_{1})_{fb} \end{array}\right] \f]
106  *
107  * where \f$\mathbf{S}_{1}\f$ is the local Schur complement matrix for
108  * each element and the subscript denotes the portion of the Schur
109  * complement associated with the vertex, edge and face blocks
110  * respectively.
111  */
113  {
114  ExpListSharedPtr expList = m_linsys.lock()->GetLocMat().lock();
116  DNekScalBlkMatSharedPtr loc_mat;
117  DNekScalMatSharedPtr bnd_mat;
118 
119  int nel, i, j, k, n, cnt, gId;
120  int meshVertId, meshEdgeId, meshFaceId;
121 
122  const int nExp = expList->GetExpSize();
123  const int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
124 
125  // Grab periodic geometry information.
127  PeriodicMap periodicVerts, periodicEdges, periodicFaces;
128  expList->GetPeriodicEntities(
129  periodicVerts, periodicEdges, periodicFaces);
130 
131  // The vectors below are of size 3 to have separate storage for
132  // vertices, edges and faces.
133 
134  // Maps from geometry ID to the matrix representing the extracted
135  // portion of S_1. For example idMats[2] folds the S_1 face blocks.
136  vector<map<int, vector<NekDouble> > > idMats(3);
137 
138  // Maps from the global ID, as obtained from AssemblyMapCG's
139  // localToGlobalMap, to the geometry ID.
140  vector<map<int, int> > gidMeshIds(3);
141 
142  // Maps from the global ID to the number of degrees of freedom for
143  // this geometry object.
144  vector<map<int, int> > gidDofs(3);
145 
146  // Array containing maximum information needed for the universal
147  // numbering later. For i = 0,1,2 for each geometry dimension:
148  //
149  // maxVertIds[2*i] = maximum geometry ID at dimension i
150  // maxVertIds[2*i+1] = maximum number of degrees of freedom for all
151  // elements of dimension i.
152  Array<OneD, int> maxVertIds(6, -1);
153 
154  // Iterator for idMats members.
155  map<int, vector<NekDouble> >::iterator gIt;
156 
157  // Figure out mapping from each elemental contribution to offset in
158  // (vert,edge,face) triples.
159  for (cnt = n = 0; n < nExp; ++n)
160  {
161  nel = expList->GetOffset_Elmt_Id(n);
162  exp = expList->GetExp(nel);
163 
164  // Grab reference to local Schur complement matrix.
165  DNekScalMatSharedPtr schurMat =
166  m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
167 
168  // Process vertices to extract relevant portion of the Schur
169  // complement matrix.
170  for (i = 0; i < exp->GetNverts(); ++i)
171  {
172  meshVertId = exp->GetGeom()->GetVid(i);
173  int locId = exp->GetVertexMap(i);
174 
175  // Get the global ID of this vertex.
176  gId = m_locToGloMap->GetLocalToGlobalMap(
177  cnt + locId) - nDirBnd;
178 
179  // Ignore all Dirichlet vertices.
180  if (gId < 0)
181  {
182  continue;
183  }
184 
185  gidDofs[0][gId] = 1;
186 
187  // Extract vertex value from Schur complement matrix.
188  NekDouble vertVal = (*schurMat)(locId,locId);
189 
190  // See if we have processed this vertex from another
191  // element.
192  gIt = idMats[0].find(gId);
193 
194  if (gIt == idMats[0].end())
195  {
196  // If not then put our 'matrix' inside idMats.
197  idMats[0][gId] = vector<NekDouble>(1, vertVal);
198  }
199  else
200  {
201  // Otherwise combine with the value that is already
202  // there (i.e. do assembly on this degree of freedom).
203  gIt->second[0] += vertVal;
204  }
205 
206  // Now check to see if the vertex is periodic. If it is,
207  // then we change meshVertId to be the minimum of all the
208  // other periodic vertices, so that we don't end up
209  // duplicating the matrix in our final block matrix.
210  pIt = periodicVerts.find(meshVertId);
211  if (pIt != periodicVerts.end())
212  {
213  for (j = 0; j < pIt->second.size(); ++j)
214  {
215  meshVertId = min(meshVertId, pIt->second[j].id);
216  }
217  }
218 
219  // Finally record the other information we need into the
220  // other matrices.
221  gidMeshIds[0][gId] = meshVertId;
222  maxVertIds[0] = max(maxVertIds[0], meshVertId);
223  maxVertIds[1] = 1;
224  }
225 
226  // Process edges. This logic is mostly the same as the previous
227  // block.
228  for (i = 0; i < exp->GetNedges(); ++i)
229  {
230  meshEdgeId = exp->GetGeom()->GetEid(i);
231 
232  Array<OneD, unsigned int> bmap, bmap2;
234  StdRegions::Orientation edgeOrient =
235  exp->GetGeom()->GetEorient(i);
236 
237  // Check if this edge is periodic. We may need to flip
238  // orientation if it is.
239  pIt = periodicEdges.find(meshEdgeId);
240  if (pIt != periodicEdges.end())
241  {
242  pair<int, StdRegions::Orientation> idOrient =
244  meshEdgeId, edgeOrient, pIt->second);
245  meshEdgeId = idOrient.first;
246  edgeOrient = idOrient.second;
247  }
248 
249  // Grab edge interior map, and the edge inverse boundary
250  // map, so that we can extract this edge from the Schur
251  // complement matrix.
252  exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
253  bmap2 = exp->GetEdgeInverseBoundaryMap(i);
254 
255  // Allocate temporary storage for the extracted edge matrix.
256  const int nEdgeCoeffs = bmap.num_elements();
257  vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
258 
259  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
260 
261  for (j = 0; j < nEdgeCoeffs; ++j)
262  {
263  // We record the minimum ID from the edge for our
264  // maps. This follows the logic that the assembly map
265  // ordering will always give us a contiguous ordering of
266  // global degrees of freedom for edge interior
267  // coefficients.
268  gId = min(gId,
269  m_locToGloMap->GetLocalToGlobalMap(
270  cnt + bmap[j])
271  - nDirBnd);
272 
273  // Ignore Dirichlet edges.
274  if (gId < 0)
275  {
276  continue;
277  }
278 
279  const NekDouble sign1 = sign[j];
280 
281  // Extract this edge, along with sign array for assembly
282  // later.
283  for (k = 0; k < nEdgeCoeffs; ++k)
284  {
285  tmpStore[k+j*nEdgeCoeffs] =
286  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
287  }
288  }
289 
290  if (gId < 0)
291  {
292  continue;
293  }
294 
295  gidDofs[1][gId] = nEdgeCoeffs;
296 
297  // Assemble this edge matrix with another one, if it exists.
298  gIt = idMats[1].find(gId);
299  if (gIt == idMats[1].end())
300  {
301  idMats[1][gId] = tmpStore;
302  }
303  else
304  {
305  ASSERTL1(tmpStore.size() == gIt->second.size(),
306  "Number of modes mismatch");
307  Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
308  &tmpStore[0], 1, &gIt->second[0], 1);
309  }
310 
311  gidMeshIds[1][gId] = meshEdgeId;
312  maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
313  maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
314  }
315 
316  // Process faces. This logic is mostly the same as the previous
317  // block.
318  for (i = 0; i < exp->GetNfaces(); ++i)
319  {
320  meshFaceId = exp->GetGeom()->GetFid(i);
321 
322  Array<OneD, unsigned int> bmap, bmap2;
324  StdRegions::Orientation faceOrient =
325  exp->GetGeom()->GetForient(i);
326 
327  // Check if this face is periodic. We may need to flip
328  // orientation if it is.
329  pIt = periodicFaces.find(meshFaceId);
330  if (pIt != periodicFaces.end())
331  {
332  meshFaceId = min(meshFaceId, pIt->second[0].id);
333  faceOrient = DeterminePeriodicFaceOrient(
334  faceOrient, pIt->second[0].orient);
335  }
336 
337  exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
338  bmap2 = exp->GetFaceInverseBoundaryMap(i);
339 
340  // Allocate temporary storage for the extracted face matrix.
341  const int nFaceCoeffs = bmap.num_elements();
342  vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
343 
344  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
345 
346  for (j = 0; j < nFaceCoeffs; ++j)
347  {
348  gId = min(gId,
349  m_locToGloMap->GetLocalToGlobalMap(
350  cnt + bmap[j])
351  - nDirBnd);
352 
353  // Ignore Dirichlet faces.
354  if (gId < 0)
355  {
356  continue;
357  }
358 
359  const NekDouble sign1 = sign[j];
360 
361  // Extract this face, along with sign array for assembly
362  // later.
363  for (k = 0; k < nFaceCoeffs; ++k)
364  {
365  tmpStore[k+j*nFaceCoeffs] =
366  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
367  }
368  }
369 
370  if (gId < 0)
371  {
372  continue;
373  }
374 
375  gidDofs[2][gId] = nFaceCoeffs;
376 
377  // Assemble this face matrix with another one, if it exists.
378  gIt = idMats[2].find(gId);
379  if (gIt == idMats[2].end())
380  {
381  idMats[2][gId] = tmpStore;
382  }
383  else
384  {
385  ASSERTL1(tmpStore.size() == gIt->second.size(),
386  "Number of modes mismatch");
387  Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
388  &tmpStore[0], 1, &gIt->second[0], 1);
389  }
390 
391  gidMeshIds[2][gId] = meshFaceId;
392  maxVertIds[4] = max(maxVertIds[4], meshFaceId);
393  maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
394  }
395 
396  cnt += exp->GetNcoeffs();
397  }
398 
399  // Perform a reduction to find maximum vertex, edge and face
400  // geometry IDs.
401  m_comm = expList->GetSession()->GetComm()->GetRowComm();
402  m_comm->AllReduce(maxVertIds, LibUtilities::ReduceMax);
403 
404  // Concatenate all matrices into contiguous storage and figure out
405  // universal ID numbering.
406  vector<NekDouble> storageBuf;
407  vector<long> globalToUniversal;
408 
409  for (i = 0, cnt = 1; i < 3; ++i)
410  {
411  const int maxDofs = maxVertIds[2*i+1];
412 
413  // Note that iterating over the map uses the property that keys
414  // are accessed in order of ascending order, putting everything
415  // in the right order for the global system.
416  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
417  {
418  // Copy matrix into storage.
419  storageBuf.insert(storageBuf.end(),
420  gIt->second.begin(), gIt->second.end());
421 
422  // Get mesh ID from global ID number.
423  ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
424  "Unable to find global ID " +
425  boost::lexical_cast<string>(gIt->first) +
426  " inside map");
427  meshVertId = gidMeshIds[i][gIt->first];
428 
429  for (j = 0; j < gIt->second.size(); ++j)
430  {
431  globalToUniversal.push_back(
432  cnt + meshVertId*maxDofs*maxDofs + j);
433  }
434 
435  // Free up the temporary storage.
436  gIt->second.clear();
437  }
438 
439  cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
440  }
441 
442  ASSERTL1(storageBuf.size() == globalToUniversal.size(),
443  "Storage buffer and global to universal map size does "
444  "not match");
445 
446  Array<OneD, NekDouble> storageData(
447  storageBuf.size(), &storageBuf[0]);
448  Array<OneD, long> globalToUniversalMap(
449  globalToUniversal.size(), &globalToUniversal[0]);
450 
451  // Use GS to assemble data between processors.
452  Gs::gs_data *tmpGs = Gs::Init(globalToUniversalMap, m_comm);
453  Gs::Gather(storageData, Gs::gs_add, tmpGs);
454 
455  // Figure out what storage we need in the block matrix.
457  1 + idMats[1].size() + idMats[2].size());
458 
459  // Vertex block is a diagonal matrix.
460  n_blks[0] = idMats[0].size();
461 
462  // Now extract number of rows in each edge and face block from the
463  // gidDofs map.
464  cnt = 1;
465  for (i = 1; i < 3; ++i)
466  {
467  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
468  {
469  ASSERTL1(gidDofs[i].count(gIt->first) > 0,
470  "Unable to find number of degrees of freedom for "
471  "global ID " + boost::lexical_cast<string>(
472  gIt->first));
473 
474  n_blks[cnt++] = gidDofs[i][gIt->first];
475  }
476  }
477 
478  // Allocate storage for the block matrix.
480  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
481 
482  // We deal with the vertex matrix separately since all vertices can
483  // be condensed into a single, block-diagonal matrix.
485  ::AllocateSharedPtr(n_blks[0], n_blks[0], 0.0, eDIAGONAL);
486 
487  // Fill the vertex matrix with the inverse of each vertex value.
488  cnt = 0;
489  for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
490  {
491  (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
492  }
493 
494  // Put the vertex matrix in the block matrix.
495  m_blkMat->SetBlock(0,0,vertMat);
496 
497  // Finally, grab the matrices from the block storage, invert them
498  // and place them in the correct position inside the block matrix.
499  int cnt2 = 1;
500  for (i = 1; i < 3; ++i)
501  {
502  for (gIt = idMats[i].begin(); gIt != idMats[i].end();
503  ++gIt, ++cnt2)
504  {
505  int nDofs = gidDofs[i][gIt->first];
506 
508  ::AllocateSharedPtr(nDofs, nDofs);
509 
510  for (j = 0; j < nDofs; ++j)
511  {
512  for (k = 0; k < nDofs; ++k)
513  {
514  (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
515  }
516  }
517 
518  cnt += nDofs*nDofs;
519 
520  tmp->Invert();
521  m_blkMat->SetBlock(cnt2, cnt2, tmp);
522  }
523  }
524  }
525 
526  /**
527  * @brief Construct a block preconditioner for the hybridized
528  * discontinuous Galerkin system.
529  *
530  * This system is built in a similar fashion to its continuous variant
531  * found in PreconditionerBlock::BlockPreconditionerCG. In this setting
532  * however, the matrix is constructed as:
533  *
534  * \f[ M^{-1} = \diag[ (\mathbf{S_{1}})_{f}^{-1} ] \f]
535  *
536  * where each matrix is the Schur complement system restricted to a
537  * single face of the trace system.
538  */
540  {
541  boost::shared_ptr<MultiRegions::ExpList>
542  expList = ((m_linsys.lock())->GetLocMat()).lock();
543  boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
545  DNekScalBlkMatSharedPtr loc_mat;
546  DNekScalMatSharedPtr bnd_mat;
547 
548  AssemblyMapDGSharedPtr asmMap = boost::dynamic_pointer_cast<
550 
551  int i, j, k, n, cnt, cnt2;
552 
553  // Figure out number of Dirichlet trace elements
554  int nTrace = expList->GetTrace()->GetExpSize();
555  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
556 
557  for (cnt = n = 0; n < nTrace; ++n)
558  {
559  if (cnt >= nDir)
560  {
561  break;
562  }
563 
564  cnt += trace->GetExp(n)->GetNcoeffs();
565  }
566 
567  nDir = n;
568 
569  // Allocate storage for block matrix. Need number of unique faces in
570  // trace space.
571  int maxTraceSize = -1;
572  map<int, int> arrayOffsets;
573 
574  for (cnt = 0, n = nDir; n < nTrace; ++n)
575  {
576  int nCoeffs = trace->GetExp(n)->GetNcoeffs();
577  int nCoeffs2 = nCoeffs * nCoeffs;
578 
579  arrayOffsets[n] = cnt;
580  cnt += nCoeffs2;
581 
582  if (maxTraceSize < nCoeffs2)
583  {
584  maxTraceSize = nCoeffs2;
585  }
586  }
587 
588  // Find maximum trace size.
589  m_comm = expList->GetSession()->GetComm()->GetRowComm();
590  m_comm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
591 
592  // Zero matrix storage.
593  Array<OneD, NekDouble> tmpStore(cnt, 0.0);
594 
595  // Assemble block matrices for each trace element.
596  for (cnt = n = 0; n < expList->GetExpSize(); ++n)
597  {
598  int elmt = expList->GetOffset_Elmt_Id(n);
599  locExpansion = expList->GetExp(elmt);
600 
602  asmMap->GetElmtToTrace()[elmt];
603 
604  // Block matrix (lambda)
605  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
606  bnd_mat = loc_mat->GetBlock(0,0);
607 
608  int nFacets = locExpansion->GetNumBases() == 2 ?
609  locExpansion->GetNedges() : locExpansion->GetNfaces();
610 
611  for (cnt2 = i = 0; i < nFacets; ++i)
612  {
613  int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
614  int elmtId = elmtToTraceMap[i]->GetElmtId ();
615 
616  // Ignore Dirichlet edges.
617  if (elmtId < nDir)
618  {
619  cnt += nCoeffs;
620  cnt2 += nCoeffs;
621  continue;
622  }
623 
624  NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
625 
626  for (j = 0; j < nCoeffs; ++j)
627  {
628  NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
629  cnt + j);
630  for (k = 0; k < nCoeffs; ++k)
631  {
632  NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
633  cnt + k);
634  off[j*nCoeffs + k] +=
635  (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
636  }
637  }
638 
639  cnt += nCoeffs;
640  cnt2 += nCoeffs;
641  }
642  }
643 
644  // Set up IDs for universal numbering.
645  Array<OneD, long> uniIds(tmpStore.num_elements());
646  for (cnt = 0, n = nDir; n < nTrace; ++n)
647  {
648  LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
649  int nCoeffs = traceExp->GetNcoeffs();
650  int geomId = traceExp->GetGeom()->GetGlobalID();
651 
652  for (i = 0; i < nCoeffs*nCoeffs; ++i)
653  {
654  uniIds[cnt++] = geomId * maxTraceSize + i + 1;
655  }
656  }
657 
658  // Assemble matrices across partitions.
659  Gs::gs_data *gsh = Gs::Init(uniIds, m_comm);
660  Gs::Gather(tmpStore, Gs::gs_add, gsh);
661 
662  // Set up diagonal block matrix
663  Array<OneD, unsigned int> n_blks(nTrace - nDir);
664  for (n = 0; n < nTrace - nDir; ++n)
665  {
666  n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
667  }
668 
670  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
671 
672  for (n = 0; n < nTrace - nDir; ++n)
673  {
674  int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
676  ::AllocateSharedPtr(nCoeffs, nCoeffs);
677  NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
678 
679  for (i = 0; i < nCoeffs; ++i)
680  {
681  for (j = 0; j < nCoeffs; ++j)
682  {
683  (*tmp)(i,j) = off[i*nCoeffs + j];
684  }
685  }
686 
687  tmp->Invert();
688  m_blkMat->SetBlock(n, n, tmp);
689  }
690  }
691 
692  /**
693  * @brief Apply preconditioner to \p pInput and store the result in \p
694  * pOutput.
695  */
697  const Array<OneD, NekDouble>& pInput,
698  Array<OneD, NekDouble>& pOutput)
699  {
700  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
701  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
702  int nNonDir = nGlobal-nDir;
703  DNekBlkMat &M = (*m_blkMat);
704  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
705  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
706  z = M * r;
707  }
708  }
709 }
boost::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:49
LibUtilities::CommSharedPtr m_comm
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Apply preconditioner to pInput and store the result in pOutput.
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:218
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::map< int, vector< PeriodicEntity > > PeriodicMap
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:150
static std::string className
Name of class.
PreconFactory & GetPreconFactory()
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
double NekDouble
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.
static PreconditionerSharedPtr create(const boost::shared_ptr< GlobalLinSys > &plinsys, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
Describe a linear system.
boost::shared_ptr< AssemblyMap > m_locToGloMap
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
const boost::weak_ptr< GlobalLinSys > m_linsys
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
PreconditionerBlock(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215