MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_VariableDofLaplacianFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Tobias Wiesner (tawiesn@sandia.gov)
42 // Ray Tuminaro (rstumin@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_
48 #define PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_
49 
50 
51 #include "MueLu_Monitor.hpp"
52 
54 
55 namespace MueLu {
56 
57  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
59  RCP<ParameterList> validParamList = rcp(new ParameterList());
60 
61  validParamList->set< double > ("Advanced Dirichlet: threshold", 1e-5, "Drop tolerance for Dirichlet detection");
62  validParamList->set< double > ("Variable DOF amalgamation: threshold", 1.8e-9, "Drop tolerance for amalgamation process");
63  validParamList->set< int > ("maxDofPerNode", 1, "Maximum number of DOFs per node");
64 
65  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
66  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for Coordinates");
67 
68  return validParamList;
69  }
70 
71  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
73 
74  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
76  Input(currentLevel, "A");
77  Input(currentLevel, "Coordinates");
78 
79  //if (currentLevel.GetLevelID() == 0) // TODO check for finest level (special treatment)
80  if (currentLevel.IsAvailable("DofPresent", NoFactory::get())) {
81  currentLevel.DeclareInput("DofPresent", NoFactory::get(), this);
82  }
83  }
84 
85  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
87  FactoryMonitor m(*this, "Build", currentLevel);
88  typedef Teuchos::ScalarTraits<SC> STS;
89 
90  const ParameterList & pL = GetParameterList();
91 
92  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
93 
94  Teuchos::RCP< const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
95  Xpetra::UnderlyingLib lib = A->getRowMap()->lib();
96 
98  RCP<dxMV> Coords = Get< RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > >(currentLevel, "Coordinates");
99 
100  int maxDofPerNode = pL.get<int>("maxDofPerNode");
101  Scalar dirDropTol = Teuchos::as<Scalar>(pL.get<double>("Advanced Dirichlet: threshold")); // "ML advnaced Dirichlet: threshold"
102  Scalar amalgDropTol = Teuchos::as<Scalar>(pL.get<double>("Variable DOF amalgamation: threshold")); //"variable DOF amalgamation: threshold")
103 
104  bool bHasZeroDiagonal = false;
106 
107  // check availability of DofPresent array
109  if (currentLevel.IsAvailable("DofPresent", NoFactory::get())) {
110  dofPresent = currentLevel.Get< Teuchos::ArrayRCP<LocalOrdinal> >("DofPresent", NoFactory::get());
111  } else {
112  // TAW: not sure about size of array. We cannot determine the expected size in the non-padded case correctly...
113  dofPresent = Teuchos::ArrayRCP<LocalOrdinal>(A->getRowMap()->getNodeNumElements(),1);
114  }
115 
116  // map[k] indicates that the kth dof in the variable dof matrix A would
117  // correspond to the map[k]th dof in the padded system. If, i.e., it is
118  // map[35] = 39 then dof no 35 in the variable dof matrix A corresponds to
119  // row map id 39 in an imaginary padded matrix Apadded.
120  // The padded system is never built but would be the associated matrix if
121  // every node had maxDofPerNode dofs.
122  std::vector<LocalOrdinal> map(A->getNodeNumRows());
123  this->buildPaddedMap(dofPresent, map, A->getNodeNumRows());
124 
125  // map of size of number of DOFs containing local node id (dof id -> node id, inclusive ghosted dofs/nodes)
126  std::vector<LocalOrdinal> myLocalNodeIds(A->getColMap()->getNodeNumElements()); // possible maximum (we need the ghost nodes, too)
127 
128  // assign the local node ids for the ghosted nodes
129  size_t nLocalNodes, nLocalPlusGhostNodes;
130  this->assignGhostLocalNodeIds(A->getRowMap(), A->getColMap(), myLocalNodeIds, map, maxDofPerNode, nLocalNodes, nLocalPlusGhostNodes, comm);
131 
132  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)," ",0,false,10,false, true);
133 
134  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dofPresent.size()) != Teuchos::as<size_t>(nLocalNodes * maxDofPerNode),MueLu::Exceptions::RuntimeError,"VariableDofLaplacianFactory: size of provided DofPresent array is " << dofPresent.size() << " but should be " << nLocalNodes * maxDofPerNode << " on the current processor.");
135 
136  // put content of assignGhostLocalNodeIds here...
137 
138  // fill nodal maps
139 
140  Teuchos::ArrayView< const GlobalOrdinal > myGids = A->getColMap()->getNodeElementList();
141 
142  // vector containing row/col gids of amalgamated matrix (with holes)
143 
144  size_t nLocalDofs = A->getRowMap()->getNodeNumElements();
145  size_t nLocalPlusGhostDofs = A->getColMap()->getNodeNumElements();
146 
147  // myLocalNodeIds (dof -> node)
148 
149  Teuchos::Array<GlobalOrdinal> amalgRowMapGIDs(nLocalNodes);
150  Teuchos::Array<GlobalOrdinal> amalgColMapGIDs(nLocalPlusGhostNodes);
151 
152  // initialize
153  size_t count = 0;
154  if (nLocalDofs > 0) {
155  amalgRowMapGIDs[count] = myGids[0];
156  amalgColMapGIDs[count] = myGids[0];
157  count++;
158  }
159 
160  for(size_t i = 1; i < nLocalDofs; i++) {
161  if (myLocalNodeIds[i] != myLocalNodeIds[i-1]) {
162  amalgRowMapGIDs[count] = myGids[i];
163  amalgColMapGIDs[count] = myGids[i];
164  count++;
165  }
166  }
167 
168  RCP<GOVector> tempAmalgColVec = GOVectorFactory::Build(A->getDomainMap());
169  Teuchos::ArrayRCP<GlobalOrdinal> tempAmalgColVecData = tempAmalgColVec->getDataNonConst(0);
170  for (size_t i = 0; i < A->getDomainMap()->getNodeNumElements(); i++)
171  tempAmalgColVecData[i] = amalgColMapGIDs[ myLocalNodeIds[i]];
172 
173  RCP<GOVector> tempAmalgColVecTarget = GOVectorFactory::Build(A->getColMap());
174  Teuchos::RCP<Import> dofImporter = ImportFactory::Build(A->getDomainMap(), A->getColMap());
175  tempAmalgColVecTarget->doImport(*tempAmalgColVec, *dofImporter, Xpetra::INSERT);
176  Teuchos::ArrayRCP<const GlobalOrdinal> tempAmalgColVecBData = tempAmalgColVecTarget->getData(0);
177 
178  // copy from dof vector to nodal vector
179  for (size_t i = 0; i < myLocalNodeIds.size(); i++)
180  amalgColMapGIDs[ myLocalNodeIds[i]] = tempAmalgColVecBData[i];
181 
182  Teuchos::RCP<Map> amalgRowMap = MapFactory::Build(lib,
184  amalgRowMapGIDs(), //View,
185  A->getRowMap()->getIndexBase(),
186  comm);
187 
188  Teuchos::RCP<Map> amalgColMap = MapFactory::Build(lib,
190  amalgColMapGIDs(), //View,
191  A->getRangeMap()->getIndexBase(),
192  comm);
193 
194  // end fill nodal maps
195 
196 
197  // start variable dof amalgamation
198 
199  Teuchos::RCP<CrsMatrixWrap> Awrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(A);
200  Teuchos::RCP<CrsMatrix> Acrs = Awrap->getCrsMatrix();
201  //Acrs->describe(*fancy, Teuchos::VERB_EXTREME);
202 
203  Teuchos::ArrayRCP<const size_t> rowptr(Acrs->getNodeNumRows());
204  Teuchos::ArrayRCP<const LocalOrdinal> colind(Acrs->getNodeNumEntries());
205  Teuchos::ArrayRCP<const Scalar> values(Acrs->getNodeNumEntries());
206  Acrs->getAllValues(rowptr, colind, values);
207 
208 
209  // create arrays for amalgamated matrix
210  Teuchos::ArrayRCP<size_t> amalgRowPtr(nLocalNodes+1);
211  Teuchos::ArrayRCP<LocalOrdinal> amalgCols(rowptr[rowptr.size()-1]);
212 
213  size_t nNonZeros = 0;
214  std::vector<bool> isNonZero(nLocalPlusGhostDofs,false);
215  std::vector<size_t> nonZeroList(nLocalPlusGhostDofs); // ???
216 
217  // also used in DetectDirichletExt
218  Teuchos::RCP<Vector> diagVecUnique = VectorFactory::Build(A->getRowMap());
219  Teuchos::RCP<Vector> diagVec = VectorFactory::Build(A->getColMap());
220  A->getLocalDiagCopy(*diagVecUnique);
221  diagVec->doImport(*diagVecUnique, *dofImporter, Xpetra::INSERT);
222  Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
223 
224  LocalOrdinal oldBlockRow = 0;
225  LocalOrdinal blockRow = 0;
226  LocalOrdinal blockColumn = 0;
227 
228  size_t newNzs = 0;
229  amalgRowPtr[0] = newNzs;
230 
231  bool doNotDrop = false;
232  if (amalgDropTol == Teuchos::ScalarTraits<Scalar>::zero()) doNotDrop = true;
233  if (values.size() == 0) doNotDrop = true;
234 
235  for(decltype(rowptr.size()) i = 0; i < rowptr.size()-1; i++) {
236  blockRow = std::floor<LocalOrdinal>( map[i] / maxDofPerNode);
237  if (blockRow != oldBlockRow) {
238  // zero out info recording nonzeros in oldBlockRow
239  for(size_t j = 0; j < nNonZeros; j++) isNonZero[nonZeroList[j]] = false;
240  nNonZeros = 0;
241  amalgRowPtr[blockRow] = newNzs; // record start of next row
242  }
243  for (size_t j = rowptr[i]; j < rowptr[i+1]; j++) {
244  if(doNotDrop == true ||
245  ( STS::magnitude(values[j] / STS::magnitude(sqrt(STS::magnitude(diagVecData[i]) * STS::magnitude(diagVecData[colind[j]]))) ) >= STS::magnitude(amalgDropTol) )) {
246  blockColumn = myLocalNodeIds[colind[j]];
247  if(isNonZero[blockColumn] == false) {
248  isNonZero[blockColumn] = true;
249  nonZeroList[nNonZeros++] = blockColumn;
250  amalgCols[newNzs++] = blockColumn;
251  }
252  }
253  }
254  oldBlockRow = blockRow;
255  }
256  amalgRowPtr[blockRow+1] = newNzs;
257 
258  TEUCHOS_TEST_FOR_EXCEPTION((blockRow+1 != Teuchos::as<LO>(nLocalNodes)) && (nLocalNodes !=0), MueLu::Exceptions::RuntimeError, "VariableDofsPerNodeAmalgamation: error, computed # block rows (" << blockRow+1 <<") != nLocalNodes (" << nLocalNodes <<")");
259 
260  amalgCols.resize(amalgRowPtr[nLocalNodes]);
261 
262  // end variableDofAmalg
263 
264  // begin rm differentDofsCrossings
265 
266  // Remove matrix entries (i,j) where the ith node and the jth node have
267  // different dofs that are 'present'
268  // Specifically, on input:
269  // dofPresent[i*maxDofPerNode+k] indicates whether or not the kth
270  // dof at the ith node is present in the
271  // variable dof matrix (e.g., the ith node
272  // has an air pressure dof). true means
273  // the dof is present while false means it
274  // is not.
275  // We create a unique id for the ith node (i.e. uniqueId[i]) via
276  // sum_{k=0 to maxDofPerNode-1} dofPresent[i*maxDofPerNode+k]*2^k
277  // and use this unique idea to remove entries (i,j) when uniqueId[i]!=uniqueId[j]
278 
279  Teuchos::ArrayRCP<LocalOrdinal> uniqueId(nLocalPlusGhostNodes); // unique id associated with DOF
280  std::vector<bool> keep(amalgRowPtr[amalgRowPtr.size()-1],true); // keep connection associated with node
281 
282  size_t ii = 0; // iteration index for present dofs
283  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
284  LocalOrdinal temp = 1; // basis for dof-id
285  uniqueId[i] = 0;
286  for (decltype(maxDofPerNode) j = 0; j < maxDofPerNode; j++) {
287  if (dofPresent[ii++]) uniqueId[i] += temp; // encode dof to be present
288  temp = temp * 2; // check next dof
289  }
290  }
291 
292  Teuchos::RCP<Import> nodeImporter = ImportFactory::Build(amalgRowMap, amalgColMap);
293 
296 
297  Teuchos::ArrayRCP< LocalOrdinal > nodeIdSrcData = nodeIdSrc->getDataNonConst(0);
298  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
299  nodeIdSrcData[i] = uniqueId[i];
300  }
301 
302  nodeIdTarget->doImport(*nodeIdSrc, *nodeImporter, Xpetra::INSERT);
303 
304  Teuchos::ArrayRCP< const LocalOrdinal > nodeIdTargetData = nodeIdTarget->getData(0);
305  for(decltype(uniqueId.size()) i = 0; i < uniqueId.size(); i++) {
306  uniqueId[i] = nodeIdTargetData[i];
307  }
308 
309  // nodal comm uniqueId, myLocalNodeIds
310 
311  // uniqueId now should contain ghosted data
312 
313  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
314  for(size_t j = amalgRowPtr[i]; j < amalgRowPtr[i+1]; j++) {
315  if (uniqueId[i] != uniqueId[amalgCols[j]]) keep [j] = false;
316  }
317  }
318 
319  // squeeze out hard-coded zeros from CSR arrays
320  Teuchos::ArrayRCP<Scalar> amalgVals;
321  this->squeezeOutNnzs(amalgRowPtr,amalgCols,amalgVals,keep);
322 
324  RCP<dxMV> ghostedCoords = dxMVf::Build(amalgColMap,Coords->getNumVectors());
325 
326  TEUCHOS_TEST_FOR_EXCEPTION(amalgRowMap->getNodeNumElements() != Coords->getMap()->getNodeNumElements(), MueLu::Exceptions::RuntimeError, "MueLu::VariableDofLaplacianFactory: the number of Coordinates and amalgamated nodes is inconsistent.");
327 
328  // Coords might live on a special nodeMap with consecutive ids (the natural numbering)
329  // The amalgRowMap might have the same number of entries, but with holes in the ids.
330  // e.g. 0,3,6,9,... as GIDs.
331  // We need the ghosted Coordinates in the buildLaplacian routine. But we access the data
332  // through getData only, i.e., the global ids are not interesting as long as we do not change
333  // the ordering of the entries
334  Coords->replaceMap(amalgRowMap);
335  ghostedCoords->doImport(*Coords, *nodeImporter, Xpetra::INSERT);
336 
337  Teuchos::ArrayRCP<Scalar> lapVals(amalgRowPtr[nLocalNodes]);
338  this->buildLaplacian(amalgRowPtr, amalgCols, lapVals, Coords->getNumVectors(), ghostedCoords);
339 
340  // sort column GIDs
341  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
342  size_t j = amalgRowPtr[i];
343  this->MueLu_az_sort<LocalOrdinal>(&(amalgCols[j]), amalgRowPtr[i+1] - j, NULL, &(lapVals[j]));
344  }
345 
346  // Caluclate status array for next level
347  Teuchos::Array<char> status(nLocalNodes * maxDofPerNode);
348 
349  // dir or not Teuchos::ArrayRCP<const bool> dirOrNot
350  for(decltype(status.size()) i = 0; i < status.size(); i++) status[i] = 's';
351  for(decltype(status.size()) i = 0; i < status.size(); i++) {
352  if(dofPresent[i] == false) status[i] = 'p';
353  }
354  if(dirOrNot.size() > 0) {
355  for(decltype(map.size()) i = 0; i < map.size(); i++) {
356  if(dirOrNot[i] == true){
357  status[map[i]] = 'd';
358  }
359  }
360  }
361  Set(currentLevel,"DofStatus",status);
362 
363  // end status array
364 
365  Teuchos::RCP<CrsMatrix> lapCrsMat = CrsMatrixFactory::Build(amalgRowMap, amalgColMap, 10); // TODO better approx for max nnz per row
366 
367  for (size_t i = 0; i < nLocalNodes; i++) {
368  lapCrsMat->insertLocalValues(i, amalgCols.view(amalgRowPtr[i],amalgRowPtr[i+1]-amalgRowPtr[i]),
369  lapVals.view(amalgRowPtr[i],amalgRowPtr[i+1]-amalgRowPtr[i]));
370  }
371  lapCrsMat->fillComplete(amalgRowMap,amalgRowMap);
372 
373  //lapCrsMat->describe(*fancy, Teuchos::VERB_EXTREME);
374 
375  Teuchos::RCP<Matrix> lapMat = Teuchos::rcp(new CrsMatrixWrap(lapCrsMat));
376  Set(currentLevel,"A",lapMat);
377  }
378 
379  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
381  TEUCHOS_TEST_FOR_EXCEPTION(numdim != 2 && numdim !=3, MueLu::Exceptions::RuntimeError,"buildLaplacian only works for 2d or 3d examples. numdim = " << numdim);
382 
383  if(numdim == 2) { // 2d
386 
387  for(decltype(rowPtr.size()) i = 0; i < rowPtr.size() - 1; i++) {
389  LocalOrdinal diag = -1;
390  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
391  if(cols[j] != Teuchos::as<LO>(i)){
392  vals[j] = std::sqrt( (x[i]-x[cols[j]]) * (x[i]-x[cols[j]]) +
393  (y[i]-y[cols[j]]) * (y[i]-y[cols[j]]) );
394  TEUCHOS_TEST_FOR_EXCEPTION(vals[j] == Teuchos::ScalarTraits<Scalar>::zero(), MueLu::Exceptions::RuntimeError, "buildLaplacian: error, " << i << " and " << cols[j] << " have same coordinates: " << x[i] << " and " << y[i]);
395  vals[j] = -Teuchos::ScalarTraits<SC>::one()/vals[j];
396  sum = sum - vals[j];
397  }
398  else diag = j;
399  }
401  TEUCHOS_TEST_FOR_EXCEPTION(diag == -1, MueLu::Exceptions::RuntimeError, "buildLaplacian: error, row " << i << " has zero diagonal!");
402 
403  vals[diag] = sum;
404  }
405  } else { // 3d
409 
410  for(decltype(rowPtr.size()) i = 0; i < rowPtr.size() - 1; i++) {
412  LocalOrdinal diag = -1;
413  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
414  if(cols[j] != Teuchos::as<LO>(i)){
415  vals[j] = std::sqrt( (x[i]-x[cols[j]]) * (x[i]-x[cols[j]]) +
416  (y[i]-y[cols[j]]) * (y[i]-y[cols[j]]) +
417  (z[i]-z[cols[j]]) * (z[i]-z[cols[j]]) );
418 
419  TEUCHOS_TEST_FOR_EXCEPTION(vals[j] == Teuchos::ScalarTraits<Scalar>::zero(), MueLu::Exceptions::RuntimeError, "buildLaplacian: error, " << i << " and " << cols[j] << " have same coordinates: " << x[i] << " and " << y[i] << " and " << z[i]);
420 
421  vals[j] = -Teuchos::ScalarTraits<SC>::one()/vals[j];
422  sum = sum - vals[j];
423  }
424  else diag = j;
425  }
427  TEUCHOS_TEST_FOR_EXCEPTION(diag == -1, MueLu::Exceptions::RuntimeError, "buildLaplacian: error, row " << i << " has zero diagonal!");
428 
429  vals[diag] = sum;
430  }
431  }
432  }
433 
434  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
436  // get rid of nonzero entries that have 0's in them and properly change
437  // the row ptr array to reflect this removal (either vals == NULL or vals != NULL)
438  // Note, the arrays are squeezed. No memory is freed.
439 
440  size_t count = 0;
441 
442  size_t nRows = rowPtr.size()-1;
443  if(vals.size() > 0) {
444  for(size_t i = 0; i < nRows; i++) {
445  size_t newStart = count;
446  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
447  if(vals[j] != Teuchos::ScalarTraits<Scalar>::zero()) {
448  cols[count ] = cols[j];
449  vals[count++] = vals[j];
450  }
451  }
452  rowPtr[i] = newStart;
453  }
454  } else {
455  for (size_t i = 0; i < nRows; i++) {
456  size_t newStart = count;
457  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
458  if (keep[j] == true) {
459  cols[count++] = cols[j];
460  }
461  }
462  rowPtr[i] = newStart;
463  }
464  }
465  rowPtr[nRows] = count;
466  }
467 
468  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
470  size_t count = 0;
471  for (decltype(dofPresent.size()) i = 0; i < dofPresent.size(); i++)
472  if(dofPresent[i] == 1) map[count++] = Teuchos::as<LocalOrdinal>(i);
473  TEUCHOS_TEST_FOR_EXCEPTION(nDofs != count, MueLu::Exceptions::RuntimeError, "VariableDofLaplacianFactory::buildPaddedMap: #dofs in dofPresent does not match the expected value (number of rows of A): " << nDofs << " vs. " << count);
474  }
475 
476  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
477  void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::assignGhostLocalNodeIds(const Teuchos::RCP<const Map> & rowDofMap, const Teuchos::RCP<const Map> & colDofMap, std::vector<LocalOrdinal> & myLocalNodeIds, const std::vector<LocalOrdinal> & dofMap, size_t maxDofPerNode, size_t& nLocalNodes, size_t& nLocalPlusGhostNodes, Teuchos::RCP< const Teuchos::Comm< int > > comm) const {
478 
479  size_t nLocalDofs = rowDofMap->getNodeNumElements();
480  size_t nLocalPlusGhostDofs = colDofMap->getNodeNumElements(); // TODO remove parameters
481 
482  // create importer for dof-based information
483  Teuchos::RCP<Import> importer = ImportFactory::Build(rowDofMap, colDofMap);
484 
485  // create a vector living on column map of A (dof based)
486  Teuchos::RCP<LOVector> localNodeIdsTemp = LOVectorFactory::Build(rowDofMap,true);
487  Teuchos::RCP<LOVector> localNodeIds = LOVectorFactory::Build(colDofMap,true);
488 
489  // fill local dofs (padded local ids)
490  Teuchos::ArrayRCP< LocalOrdinal > localNodeIdsTempData = localNodeIdsTemp->getDataNonConst(0);
491  for(size_t i = 0; i < localNodeIdsTemp->getLocalLength(); i++)
492  localNodeIdsTempData[i] = std::floor<LocalOrdinal>( dofMap[i] / maxDofPerNode );
493 
494  localNodeIds->doImport(*localNodeIdsTemp, *importer, Xpetra::INSERT);
495  Teuchos::ArrayRCP< const LocalOrdinal > localNodeIdsData = localNodeIds->getData(0);
496 
497  // Note: localNodeIds contains local ids for the padded version as vector values
498 
499 
500  // we use Scalar instead of int as type
501  Teuchos::RCP<LOVector> myProcTemp = LOVectorFactory::Build(rowDofMap,true);
502  Teuchos::RCP<LOVector> myProc = LOVectorFactory::Build(colDofMap,true);
503 
504  // fill local dofs (padded local ids)
505  Teuchos::ArrayRCP< LocalOrdinal > myProcTempData = myProcTemp->getDataNonConst(0);
506  for(size_t i = 0; i < myProcTemp->getLocalLength(); i++)
507  myProcTempData[i] = Teuchos::as<LocalOrdinal>(comm->getRank());
508  myProc->doImport(*myProcTemp, *importer, Xpetra::INSERT);
509  Teuchos::ArrayRCP<LocalOrdinal> myProcData = myProc->getDataNonConst(0); // we have to modify the data (therefore the non-const version)
510 
511  // At this point, the ghost part of localNodeIds corresponds to the local ids
512  // associated with the current owning processor. We want to convert these to
513  // local ids associated with the processor on which these are ghosts.
514  // Thus we have to re-number them. In doing this re-numbering we must make sure
515  // that we find all ghosts with the same id & proc and assign a unique local
516  // id to this group (id&proc). To do this find, we sort all ghost entries in
517  // localNodeIds that are owned by the same processor. Then we can look for
518  // duplicates (i.e., several ghost entries corresponding to dofs with the same
519  // node id) easily and make sure these are all assigned to the same local id.
520  // To do the sorting we'll make a temporary copy of the ghosts via tempId and
521  // tempProc and sort this multiple times for each group owned by the same proc.
522 
523 
524  std::vector<size_t> location(nLocalPlusGhostDofs - nLocalDofs + 1);
525  std::vector<size_t> tempId (nLocalPlusGhostDofs - nLocalDofs + 1);
526  std::vector<size_t> tempProc(nLocalPlusGhostDofs - nLocalDofs + 1);
527 
528  size_t notProcessed = nLocalDofs; // iteration index over all ghosted dofs
529  size_t tempIndex = 0;
530  size_t first = tempIndex;
531  LocalOrdinal neighbor;
532 
533  while (notProcessed < nLocalPlusGhostDofs) {
534  neighbor = myProcData[notProcessed]; // get processor id of not-processed element
535  first = tempIndex;
536  location[tempIndex] = notProcessed;
537  tempId[tempIndex++] = localNodeIdsData[notProcessed];
538  myProcData[notProcessed] = -1 - neighbor;
539 
540  for(size_t i = notProcessed + 1; i < nLocalPlusGhostDofs; i++) {
541  if(myProcData[i] == neighbor) {
542  location[tempIndex] = i;
543  tempId[tempIndex++] = localNodeIdsData[i];
544  myProcData[i] = -1; // mark as visited
545  }
546  }
547  this->MueLu_az_sort<size_t>(&(tempId[first]), tempIndex - first, &(location[first]), NULL);
548  for(size_t i = first; i < tempIndex; i++) tempProc[i] = neighbor;
549 
550  // increment index. Find next notProcessed dof index corresponding to first non-visited element
551  notProcessed++;
552  while ( (notProcessed < nLocalPlusGhostDofs) && (myProcData[notProcessed] < 0))
553  notProcessed++;
554  }
555  TEUCHOS_TEST_FOR_EXCEPTION(tempIndex != nLocalPlusGhostDofs-nLocalDofs, MueLu::Exceptions::RuntimeError,"Number of nonzero ghosts is inconsistent.");
556 
557  // Now assign ids to all ghost nodes (giving the same id to those with the
558  // same myProc[] and the same local id on the proc that actually owns the
559  // variable associated with the ghost
560 
561  nLocalNodes = 0; // initialize return value
562  if(nLocalDofs > 0) nLocalNodes = localNodeIdsData[nLocalDofs-1] + 1;
563 
564  nLocalPlusGhostNodes = nLocalNodes; // initialize return value
565  if(nLocalDofs < nLocalPlusGhostDofs) nLocalPlusGhostNodes++; // 1st ghost node is unique (not accounted for). number will be increased later, if there are more ghost nodes
566 
567  // check if two adjacent ghost dofs correspond to different nodes. To do this,
568  // check if they are from different processors or whether they have different
569  // local node ids
570 
571  // loop over all (remaining) ghost dofs
572  for (size_t i = nLocalDofs+1; i < nLocalPlusGhostDofs; i++) {
573  size_t lagged = nLocalPlusGhostNodes-1;
574 
575  // i is a new unique ghost node (not already accounted for)
576  if ((tempId[i-nLocalDofs] != tempId[i-1-nLocalDofs]) ||
577  (tempProc[i-nLocalDofs] != tempProc[i-1-nLocalDofs]))
578  nLocalPlusGhostNodes++; // update number of ghost nodes
579  tempId[i-1-nLocalDofs] = lagged;
580  }
581  if (nLocalPlusGhostDofs > nLocalDofs)
582  tempId[nLocalPlusGhostDofs-1-nLocalDofs] = nLocalPlusGhostNodes - 1;
583 
584  // fill myLocalNodeIds array. Start with local part (not ghosted)
585  for(size_t i = 0; i < nLocalDofs; i++)
586  myLocalNodeIds[i] = std::floor<LocalOrdinal>( dofMap[i] / maxDofPerNode );
587 
588  // copy ghosted nodal ids into myLocalNodeIds
589  for(size_t i = nLocalDofs; i < nLocalPlusGhostDofs; i++)
590  myLocalNodeIds[location[i-nLocalDofs]] = tempId[i-nLocalDofs];
591 
592  }
593 
594 } /* MueLu */
595 
596 
597 #endif /* PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_ */
void buildPaddedMap(const Teuchos::ArrayRCP< const LocalOrdinal > &dofPresent, std::vector< LocalOrdinal > &map, size_t nDofs) const
void buildLaplacian(const Teuchos::ArrayRCP< size_t > &rowPtr, const Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const size_t &numdim, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > &ghostedCoords) const
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LocalOrdinal LO
size_type size() const
void DeclareInput(Level &currentLevel) const
Input.
void assignGhostLocalNodeIds(const Teuchos::RCP< const Map > &rowDofMap, const Teuchos::RCP< const Map > &colDofMap, std::vector< LocalOrdinal > &myLocalNodeIds, const std::vector< LocalOrdinal > &dofMap, size_t maxDofPerNode, size_t &nLocalNodes, size_t &nLocalPlusGhostNodes, Teuchos::RCP< const Teuchos::Comm< int > > comm) const
static const NoFactory * get()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &currentLevel) const
Build an object with this factory.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
size_type size() const
void squeezeOutNnzs(Teuchos::ArrayRCP< size_t > &rowPtr, Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const std::vector< bool > &keep) const
Node NO
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
ArrayView< T > view(size_type lowerOffset, size_type size) const