Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GlobalLinSysStaticCond.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: GlobalLinSysStaticCond.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: Implementation to linear solver using single-
33 // or multi-level static condensation
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <MultiRegions/GlobalLinSysStaticCond.h>
38 #include <LibUtilities/BasicUtils/Timer.h>
39 #include <LibUtilities/BasicUtils/ErrorUtil.hpp>
40 #include <LibUtilities/LinearAlgebra/StorageSmvBsr.hpp>
41 #include <LibUtilities/LinearAlgebra/SparseDiagBlkMatrix.hpp>
42 #include <LibUtilities/LinearAlgebra/SparseUtils.hpp>
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
48  /**
49  * @class GlobalLinSysStaticCond
50  *
51  * Solves a linear system using single- or multi-level static
52  * condensation.
53  */
54 
55  /**
56  * For a matrix system of the form @f[
57  * \left[ \begin{array}{cc}
58  * \boldsymbol{A} & \boldsymbol{B}\\
59  * \boldsymbol{C} & \boldsymbol{D}
60  * \end{array} \right]
61  * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2}
62  * \end{array}\right]
63  * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2}
64  * \end{array}\right],
65  * @f]
66  * where @f$\boldsymbol{D}@f$ and
67  * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble
68  * a static condensation system, according to a given local to global
69  * mapping. #m_linSys is constructed by AssembleSchurComplement().
70  * @param mKey Associated matrix key.
71  * @param pLocMatSys LocalMatrixSystem
72  * @param locToGloMap Local to global mapping.
73  */
74  GlobalLinSysStaticCond::GlobalLinSysStaticCond(
75  const GlobalLinSysKey &pKey,
76  const boost::weak_ptr<ExpList> &pExpList,
77  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
78  : GlobalLinSys(pKey, pExpList, pLocToGloMap),
79  m_locToGloMap (pLocToGloMap)
80  {
81  }
82 
83  void GlobalLinSysStaticCond::v_InitObject()
84  {
85  // Allocate memory for top-level structure
86  SetupTopLevel(m_locToGloMap);
87 
88  // Construct this level
89  Initialise(m_locToGloMap);
90  }
91 
92  /**
93  *
94  */
95  GlobalLinSysStaticCond::~GlobalLinSysStaticCond()
96  {
97 
98  }
99 
100 
101  /**
102  *
103  */
104  void GlobalLinSysStaticCond::v_Solve(
105  const Array<OneD, const NekDouble> &in,
106  Array<OneD, NekDouble> &out,
107  const AssemblyMapSharedPtr &pLocToGloMap,
108  const Array<OneD, const NekDouble> &dirForcing)
109  {
110  bool dirForcCalculated = (bool) dirForcing.num_elements();
111  bool atLastLevel = pLocToGloMap->AtLastLevel();
112  int scLevel = pLocToGloMap->GetStaticCondLevel();
113 
114  int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
115  int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
116  int nDirBndDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
117  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
118  int nLocBndDofs = pLocToGloMap->GetNumLocalBndCoeffs();
119  int nIntDofs = pLocToGloMap->GetNumGlobalCoeffs()
120  - nGlobBndDofs;
121 
122  Array<OneD, NekDouble> F = m_wsp + 2*nLocBndDofs;
123  Array<OneD, NekDouble> tmp;
124  if(nDirBndDofs && dirForcCalculated)
125  {
126  Vmath::Vsub(nGlobDofs,in.get(),1,dirForcing.get(),1,F.get(),1);
127  }
128  else
129  {
130  Vmath::Vcopy(nGlobDofs,in.get(),1,F.get(),1);
131  }
132 
133  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,tmp=F+nDirBndDofs,
134  eWrapper);
135  NekVector<NekDouble> F_GlobBnd(nGlobBndDofs,F,eWrapper);
136  NekVector<NekDouble> F_LocBnd(nLocBndDofs,0.0);
137  NekVector<NekDouble> F_Int(nIntDofs,tmp=F+nGlobBndDofs,eWrapper);
138 
139  NekVector<NekDouble> V_GlobBnd(nGlobBndDofs,out,eWrapper);
140  NekVector<NekDouble> V_GlobHomBnd(nGlobHomBndDofs,
141  tmp=out+nDirBndDofs,
142  eWrapper);
143  NekVector<NekDouble> V_Int(nIntDofs,tmp=out+nGlobBndDofs,eWrapper);
144  NekVector<NekDouble> V_LocBnd(nLocBndDofs,m_wsp,eWrapper);
145 
146  NekVector<NekDouble> V_GlobHomBndTmp(nGlobHomBndDofs,0.0);
147 
148  // set up normalisation factor for right hand side on first SC level
149  DNekScalBlkMatSharedPtr sc = v_PreSolve(scLevel, F_GlobBnd);
150 
151  if(nGlobHomBndDofs)
152  {
153  // construct boundary forcing
154  if( nIntDofs && ((!dirForcCalculated) && (atLastLevel)) )
155  {
156  DNekScalBlkMat &BinvD = *m_BinvD;
157  DNekScalBlkMat &SchurCompl = *sc;
158 
159  // include dirichlet boundary forcing
160  pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
161  V_LocBnd = BinvD*F_Int + SchurCompl*V_LocBnd;
162  }
163  else if((!dirForcCalculated) && (atLastLevel))
164  {
165  // include dirichlet boundary forcing
166  DNekScalBlkMat &SchurCompl = *sc;
167  pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
168  V_LocBnd = SchurCompl*V_LocBnd;
169  }
170  else
171  {
172  DNekScalBlkMat &BinvD = *m_BinvD;
173  V_LocBnd = BinvD*F_Int;
174  }
175 
176  pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
177  nDirBndDofs);
178  F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
179 
180  // Transform from original basis to low energy
181  v_BasisTransform(F, nDirBndDofs);
182 
183  // For parallel multi-level static condensation some
184  // processors may have different levels to others. This
185  // routine receives contributions to partition vertices from
186  // those lower levels, whilst not sending anything to the
187  // other partitions, and includes them in the modified right
188  // hand side vector.
189  int lcLevel = pLocToGloMap->GetLowestStaticCondLevel();
190  if(atLastLevel && scLevel < lcLevel)
191  {
192  // If this level is not the lowest level across all
193  // processes, we must do dummy communication for the
194  // remaining levels
195  Array<OneD, NekDouble> tmp(nGlobBndDofs);
196  for (int i = scLevel; i < lcLevel; ++i)
197  {
198  Vmath::Fill(nGlobBndDofs, 0.0, tmp, 1);
199  pLocToGloMap->UniversalAssembleBnd(tmp);
200  Vmath::Vcopy(nGlobHomBndDofs,
201  tmp.get()+nDirBndDofs, 1,
202  V_GlobHomBndTmp.GetPtr().get(), 1);
203  F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
204  }
205  }
206 
207  // solve boundary system
208  if(atLastLevel)
209  {
210  Array<OneD, NekDouble> pert(nGlobBndDofs,0.0);
211  NekVector<NekDouble> Pert(nGlobBndDofs,pert,eWrapper);
212 
213  // Solve for difference from initial solution given inout;
214  SolveLinearSystem(
215  nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
216 
217  // Transform back to original basis
218  v_BasisInvTransform(pert);
219 
220  // Add back initial conditions onto difference
221  Vmath::Vadd(nGlobHomBndDofs,&out[nDirBndDofs],1,
222  &pert[nDirBndDofs],1,&out[nDirBndDofs],1);
223  }
224  else
225  {
226  m_recursiveSchurCompl->Solve(F,
227  V_GlobBnd.GetPtr(),
228  pLocToGloMap->GetNextLevelLocalToGlobalMap());
229  }
230  }
231 
232  // solve interior system
233  if(nIntDofs)
234  {
235  DNekScalBlkMat &invD = *m_invD;
236 
237  if(nGlobHomBndDofs || nDirBndDofs)
238  {
239  DNekScalBlkMat &C = *m_C;
240 
241  if(dirForcCalculated && nDirBndDofs)
242  {
243  pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd,
244  nDirBndDofs);
245  }
246  else
247  {
248  pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
249  }
250  F_Int = F_Int - C*V_LocBnd;
251  }
252 
253  V_Int = invD*F_Int;
254  }
255  }
256 
257 
258  /**
259  * If at the last level of recursion (or the only level in the case of
260  * single-level static condensation), assemble the Schur complement.
261  * For other levels, in the case of multi-level static condensation,
262  * the next level of the condensed system is computed.
263  * @param pLocToGloMap Local to global mapping.
264  */
265  void GlobalLinSysStaticCond::v_Initialise(
266  const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
267  {
268  int nLocalBnd = m_locToGloMap->GetNumLocalBndCoeffs();
269  int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
270  m_wsp = Array<OneD, NekDouble>(2*nLocalBnd + nGlobal, 0.0);
271 
272  if (pLocToGloMap->AtLastLevel())
273  {
274  v_AssembleSchurComplement(pLocToGloMap);
275  }
276  else
277  {
278  ConstructNextLevelCondensedSystem(
279  pLocToGloMap->GetNextLevelLocalToGlobalMap());
280  }
281  }
282 
283  int GlobalLinSysStaticCond::v_GetNumBlocks()
284  {
285  return m_schurCompl->GetNumberOfBlockRows();
286  }
287 
288  /**
289  * For the first level in multi-level static condensation, or the only
290  * level in the case of single-level static condensation, allocate the
291  * condensed matrices and populate them with the local matrices
292  * retrieved from the expansion list.
293  * @param
294  */
295  void GlobalLinSysStaticCond::SetupTopLevel(
296  const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
297  {
298  int n;
299  int n_exp = m_expList.lock()->GetNumElmts();
300 
301  const Array<OneD,const unsigned int>& nbdry_size
302  = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
303  const Array<OneD,const unsigned int>& nint_size
304  = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
305 
306  // Setup Block Matrix systems
307  MatrixStorage blkmatStorage = eDIAGONAL;
308  m_schurCompl = MemoryManager<DNekScalBlkMat>
309  ::AllocateSharedPtr(nbdry_size, nbdry_size, blkmatStorage);
310  m_BinvD = MemoryManager<DNekScalBlkMat>
311  ::AllocateSharedPtr(nbdry_size, nint_size , blkmatStorage);
312  m_C = MemoryManager<DNekScalBlkMat>
313  ::AllocateSharedPtr(nint_size , nbdry_size, blkmatStorage);
314  m_invD = MemoryManager<DNekScalBlkMat>
315  ::AllocateSharedPtr(nint_size , nint_size , blkmatStorage);
316 
317  for(n = 0; n < n_exp; ++n)
318  {
319  if (m_linSysKey.GetMatrixType() ==
320  StdRegions::eHybridDGHelmBndLam)
321  {
322  DNekScalMatSharedPtr loc_mat
323  = GlobalLinSys::v_GetBlock(
324  m_expList.lock()->GetOffset_Elmt_Id(n));
325  m_schurCompl->SetBlock(n,n,loc_mat);
326  }
327  else
328  {
329  DNekScalBlkMatSharedPtr loc_schur
330  = GlobalLinSys::v_GetStaticCondBlock(
331  m_expList.lock()->GetOffset_Elmt_Id(n));
332  DNekScalMatSharedPtr t;
333  m_schurCompl->SetBlock(n, n, t = loc_schur->GetBlock(0,0));
334  m_BinvD ->SetBlock(n, n, t = loc_schur->GetBlock(0,1));
335  m_C ->SetBlock(n, n, t = loc_schur->GetBlock(1,0));
336  m_invD ->SetBlock(n, n, t = loc_schur->GetBlock(1,1));
337  }
338  }
339  }
340 
341  /**
342  *
343  */
344  void GlobalLinSysStaticCond::ConstructNextLevelCondensedSystem(
345  const AssemblyMapSharedPtr& pLocToGloMap)
346  {
347  int i,j,n,cnt;
348  DNekScalBlkMatSharedPtr blkMatrices[4];
349 
350  // Create temporary matrices within an inner-local scope to ensure
351  // any references to the intermediate storage is lost before
352  // the recursive step, rather than at the end of the routine.
353  // This allows the schur complement matrix from this level to be
354  // disposed of in the next level after use without being retained
355  // due to lingering shared pointers.
356  {
357 
358  const Array<OneD,const unsigned int>& nBndDofsPerPatch
359  = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
360  const Array<OneD,const unsigned int>& nIntDofsPerPatch
361  = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
362 
363  // STEP 1:
364  // Based upon the schur complement of the the current level we
365  // will substructure this matrix in the form
366  // -- --
367  // | A B |
368  // | C D |
369  // -- --
370  // All matrices A,B,C and D are (diagonal) blockmatrices.
371  // However, as a start we will use an array of DNekMatrices as
372  // it is too hard to change the individual entries of a
373  // DNekScalBlkMatSharedPtr.
374 
375  // In addition, we will also try to ensure that the memory of
376  // the blockmatrices will be contiguous. This will probably
377  // enhance the efficiency
378  // - Calculate the total number of entries in the blockmatrices
379  int nPatches = pLocToGloMap->GetNumPatches();
380  int nEntriesA = 0; int nEntriesB = 0;
381  int nEntriesC = 0; int nEntriesD = 0;
382 
383  for(i = 0; i < nPatches; i++)
384  {
385  nEntriesA += nBndDofsPerPatch[i]*nBndDofsPerPatch[i];
386  nEntriesB += nBndDofsPerPatch[i]*nIntDofsPerPatch[i];
387  nEntriesC += nIntDofsPerPatch[i]*nBndDofsPerPatch[i];
388  nEntriesD += nIntDofsPerPatch[i]*nIntDofsPerPatch[i];
389  }
390 
391  // Now create the DNekMatrices and link them to the memory
392  // allocated above
393  Array<OneD, DNekMatSharedPtr> substructuredMat[4]
394  = {Array<OneD, DNekMatSharedPtr>(nPatches), //Matrix A
395  Array<OneD, DNekMatSharedPtr>(nPatches), //Matrix B
396  Array<OneD, DNekMatSharedPtr>(nPatches), //Matrix C
397  Array<OneD, DNekMatSharedPtr>(nPatches)}; //Matrix D
398 
399  // Initialise storage for the matrices. We do this separately
400  // for each matrix so the matrices may be independently
401  // deallocated when no longer required.
402  Array<OneD, NekDouble> storageA(nEntriesA,0.0);
403  Array<OneD, NekDouble> storageB(nEntriesB,0.0);
404  Array<OneD, NekDouble> storageC(nEntriesC,0.0);
405  Array<OneD, NekDouble> storageD(nEntriesD,0.0);
406 
407  Array<OneD, NekDouble> tmparray;
408  PointerWrapper wType = eWrapper;
409  int cntA = 0;
410  int cntB = 0;
411  int cntC = 0;
412  int cntD = 0;
413 
414  for(i = 0; i < nPatches; i++)
415  {
416  // Matrix A
417  tmparray = storageA+cntA;
418  substructuredMat[0][i] = MemoryManager<DNekMat>
419  ::AllocateSharedPtr(nBndDofsPerPatch[i],
420  nBndDofsPerPatch[i],
421  tmparray, wType);
422  // Matrix B
423  tmparray = storageB+cntB;
424  substructuredMat[1][i] = MemoryManager<DNekMat>
425  ::AllocateSharedPtr(nBndDofsPerPatch[i],
426  nIntDofsPerPatch[i],
427  tmparray, wType);
428  // Matrix C
429  tmparray = storageC+cntC;
430  substructuredMat[2][i] = MemoryManager<DNekMat>
431  ::AllocateSharedPtr(nIntDofsPerPatch[i],
432  nBndDofsPerPatch[i],
433  tmparray, wType);
434  // Matrix D
435  tmparray = storageD+cntD;
436  substructuredMat[3][i] = MemoryManager<DNekMat>
437  ::AllocateSharedPtr(nIntDofsPerPatch[i],
438  nIntDofsPerPatch[i],
439  tmparray, wType);
440 
441  cntA += nBndDofsPerPatch[i] * nBndDofsPerPatch[i];
442  cntB += nBndDofsPerPatch[i] * nIntDofsPerPatch[i];
443  cntC += nIntDofsPerPatch[i] * nBndDofsPerPatch[i];
444  cntD += nIntDofsPerPatch[i] * nIntDofsPerPatch[i];
445  }
446 
447  // Then, project SchurComplement onto
448  // the substructured matrices of the next level
449  DNekScalBlkMatSharedPtr SchurCompl = m_schurCompl;
450  DNekScalMatSharedPtr schurComplSubMat;
451  int schurComplSubMatnRows;
452  Array<OneD, const int> patchId, dofId;
453  Array<OneD, const unsigned int> isBndDof;
454  Array<OneD, const NekDouble> sign;
455  NekDouble scale;
456 
457  for(n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
458  {
459  schurComplSubMat = SchurCompl->GetBlock(n,n);
460  schurComplSubMatnRows = schurComplSubMat->GetRows();
461 
462  scale = SchurCompl->GetBlock(n,n)->Scale();
463  Array<OneD, NekDouble> schurSubMat
464  = SchurCompl->GetBlock(n,n)->GetOwnedMatrix()->GetPtr();
465 
466  patchId = pLocToGloMap->GetPatchMapFromPrevLevel()
467  ->GetPatchId() + cnt;
468  dofId = pLocToGloMap->GetPatchMapFromPrevLevel()
469  ->GetDofId() + cnt;
470  isBndDof = pLocToGloMap->GetPatchMapFromPrevLevel()
471  ->IsBndDof() + cnt;
472  sign = pLocToGloMap->GetPatchMapFromPrevLevel()
473  ->GetSign() + cnt;
474 
475  // Set up Matrix;
476  for(i = 0; i < schurComplSubMatnRows; ++i)
477  {
478  int pId = patchId[i];
479  Array<OneD, NekDouble> subMat0
480  = substructuredMat[0][pId]->GetPtr();
481  Array<OneD, NekDouble> subMat1
482  = substructuredMat[1][patchId[i]]->GetPtr();
483  Array<OneD, NekDouble> subMat2
484  = substructuredMat[2][patchId[i]]->GetPtr();
485  Array<OneD, NekDouble> subMat3
486  = substructuredMat[3][patchId[i]]->GetPtr();
487  int subMat0rows = substructuredMat[0][pId]->GetRows();
488  int subMat1rows = substructuredMat[1][pId]->GetRows();
489  int subMat2rows = substructuredMat[2][pId]->GetRows();
490  int subMat3rows = substructuredMat[3][pId]->GetRows();
491 
492  if(isBndDof[i])
493  {
494  for(j = 0; j < schurComplSubMatnRows; ++j)
495  {
496  ASSERTL0(patchId[i]==patchId[j],
497  "These values should be equal");
498 
499  if(isBndDof[j])
500  {
501  subMat0[dofId[i]+dofId[j]*subMat0rows] +=
502  sign[i]*sign[j]*(
503  scale*schurSubMat[
504  i+j*schurComplSubMatnRows]);
505  }
506  else
507  {
508  subMat1[dofId[i]+dofId[j]*subMat1rows] +=
509  sign[i]*sign[j]*(
510  scale*schurSubMat[
511  i+j*schurComplSubMatnRows]);
512  }
513  }
514  }
515  else
516  {
517  for(j = 0; j < schurComplSubMatnRows; ++j)
518  {
519  ASSERTL0(patchId[i]==patchId[j],
520  "These values should be equal");
521 
522  if(isBndDof[j])
523  {
524  subMat2[dofId[i]+dofId[j]*subMat2rows] +=
525  sign[i]*sign[j]*(
526  scale*schurSubMat[
527  i+j*schurComplSubMatnRows]);
528  }
529  else
530  {
531  subMat3[dofId[i]+dofId[j]*subMat3rows] +=
532  sign[i]*sign[j]*(
533  scale*schurSubMat[
534  i+j*schurComplSubMatnRows]);
535  }
536  }
537  }
538  }
539  cnt += schurComplSubMatnRows;
540  }
541 
542  // STEP 2: condense the system
543  // This can be done elementally (i.e. patch per patch)
544  for(i = 0; i < nPatches; i++)
545  {
546  if(nIntDofsPerPatch[i])
547  {
548  Array<OneD, NekDouble> subMat0
549  = substructuredMat[0][i]->GetPtr();
550  Array<OneD, NekDouble> subMat1
551  = substructuredMat[1][i]->GetPtr();
552  Array<OneD, NekDouble> subMat2
553  = substructuredMat[2][i]->GetPtr();
554  int subMat0rows = substructuredMat[0][i]->GetRows();
555  int subMat1rows = substructuredMat[1][i]->GetRows();
556  int subMat2rows = substructuredMat[2][i]->GetRows();
557  int subMat2cols = substructuredMat[2][i]->GetColumns();
558 
559  // 1. D -> InvD
560  substructuredMat[3][i]->Invert();
561  // 2. B -> BInvD
562  (*substructuredMat[1][i]) = (*substructuredMat[1][i])*
563  (*substructuredMat[3][i]);
564  // 3. A -> A - BInvD*C (= schurcomplement)
565  // (*substructuredMat[0][i]) =(*substructuredMat[0][i])-
566  // (*substructuredMat[1][i])*(*substructuredMat[2][i]);
567  // Note: faster to use blas directly
568  Blas::Dgemm('N','N', subMat1rows, subMat2cols,
569  subMat2rows, -1.0, &subMat1[0], subMat1rows,
570  &subMat2[0], subMat2rows, 1.0, &subMat0[0],
571  subMat0rows);
572  }
573  }
574 
575  // STEP 3: fill the block matrices. However, do note that we
576  // first have to convert them to a DNekScalMat in order to be
577  // compatible with the first level of static condensation
578  const Array<OneD,const unsigned int>& nbdry_size
579  = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
580  const Array<OneD,const unsigned int>& nint_size
581  = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
582  MatrixStorage blkmatStorage = eDIAGONAL;
583 
584  blkMatrices[0] = MemoryManager<DNekScalBlkMat>
585  ::AllocateSharedPtr(nbdry_size, nbdry_size, blkmatStorage);
586  blkMatrices[1] = MemoryManager<DNekScalBlkMat>
587  ::AllocateSharedPtr(nbdry_size, nint_size , blkmatStorage);
588  blkMatrices[2] = MemoryManager<DNekScalBlkMat>
589  ::AllocateSharedPtr(nint_size , nbdry_size, blkmatStorage);
590  blkMatrices[3] = MemoryManager<DNekScalBlkMat>
591  ::AllocateSharedPtr(nint_size , nint_size , blkmatStorage);
592 
593  DNekScalMatSharedPtr tmpscalmat;
594  for(i = 0; i < nPatches; i++)
595  {
596  for(j = 0; j < 4; j++)
597  {
598  tmpscalmat= MemoryManager<DNekScalMat>
599  ::AllocateSharedPtr(1.0,substructuredMat[j][i]);
600  blkMatrices[j]->SetBlock(i,i,tmpscalmat);
601  }
602  }
603  }
604 
605  // We've finished with the Schur complement matrix passed to this
606  // level, so return the memory to the system. The Schur complement
607  // matrix need only be retained at the last level. Save the other
608  // matrices at this level though.
609  m_schurCompl.reset();
610 
611  m_recursiveSchurCompl = v_Recurse(
612  m_linSysKey, m_expList, blkMatrices[0], blkMatrices[1],
613  blkMatrices[2], blkMatrices[3], pLocToGloMap);
614  }
615  }
616 }
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2