MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_CoalesceDropFactory_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 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_COALESCEDROPFACTORY_DEF_HPP
47 #define MUELU_COALESCEDROPFACTORY_DEF_HPP
48 
49 #include <Xpetra_CrsGraphFactory.hpp>
50 #include <Xpetra_CrsGraph.hpp>
51 #include <Xpetra_ImportFactory.hpp>
52 #include <Xpetra_MapFactory.hpp>
53 #include <Xpetra_Map.hpp>
54 #include <Xpetra_Matrix.hpp>
55 #include <Xpetra_MultiVectorFactory.hpp>
56 #include <Xpetra_MultiVector.hpp>
57 #include <Xpetra_StridedMap.hpp>
58 #include <Xpetra_VectorFactory.hpp>
59 #include <Xpetra_Vector.hpp>
60 
61 #include "MueLu_CoalesceDropFactory_decl.hpp"
62 
63 #include "MueLu_AmalgamationFactory.hpp"
64 #include "MueLu_AmalgamationInfo.hpp"
65 #include "MueLu_Exceptions.hpp"
66 #include "MueLu_GraphBase.hpp"
67 #include "MueLu_Graph.hpp"
68 #include "MueLu_Level.hpp"
69 #include "MueLu_LWGraph.hpp"
70 #include "MueLu_MasterList.hpp"
71 #include "MueLu_Monitor.hpp"
72 #include "MueLu_PreDropFunctionBaseClass.hpp"
73 #include "MueLu_PreDropFunctionConstVal.hpp"
74 #include "MueLu_Utilities.hpp"
75 
76 namespace MueLu {
77 
78  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
79  RCP<const ParameterList> CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
80  RCP<ParameterList> validParamList = rcp(new ParameterList());
81 
82 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
83  SET_VALID_ENTRY("aggregation: drop tol");
84  SET_VALID_ENTRY("aggregation: Dirichlet threshold");
85  SET_VALID_ENTRY("aggregation: drop scheme");
86  {
87  typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
88  validParamList->getEntry("aggregation: drop scheme").setValidator(
89  rcp(new validatorType(Teuchos::tuple<std::string>("classical", "distance laplacian"), "aggregation: drop scheme")));
90  }
91 #undef SET_VALID_ENTRY
92  validParamList->set< bool > ("lightweight wrap", true, "Experimental option for lightweight graph access");
93 
94  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
95  validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory for UnAmalgamationInfo");
96  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for Coordinates");
97 
98  return validParamList;
99  }
100 
101  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
102  CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CoalesceDropFactory() : predrop_(Teuchos::null) { }
103 
104  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
105  void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &currentLevel) const {
106  Input(currentLevel, "A");
107  Input(currentLevel, "UnAmalgamationInfo");
108 
109  const ParameterList& pL = GetParameterList();
110  if (pL.get<bool>("lightweight wrap") == true) {
111  if (pL.get<std::string>("aggregation: drop scheme") == "distance laplacian")
112  Input(currentLevel, "Coordinates");
113 
114  }
115  }
116 
117  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
118  void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const {
119  FactoryMonitor m(*this, "Build", currentLevel);
120 
121  typedef Teuchos::ScalarTraits<SC> STS;
122  typedef typename STS::magnitudeType real_type;
123  typedef Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
124 
125  if (predrop_ != Teuchos::null)
126  GetOStream(Parameters0) << predrop_->description();
127 
128  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
129  RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel, "UnAmalgamationInfo");
130 
131  const ParameterList & pL = GetParameterList();
132  bool doExperimentalWrap = pL.get<bool>("lightweight wrap");
133 
134  GetOStream(Parameters0) << "lightweight wrap = " << doExperimentalWrap << std::endl;
135 
136  // decide wether to use the fast-track code path for standard maps or the somewhat slower
137  // code path for non-standard maps
138  /*bool bNonStandardMaps = false;
139  if (A->IsView("stridedMaps") == true) {
140  Teuchos::RCP<const Map> myMap = A->getRowMap("stridedMaps");
141  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
142  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
143  if (strMap->getStridedBlockId() != -1 || strMap->getOffset() > 0)
144  bNonStandardMaps = true;
145  }*/
146 
147  if (doExperimentalWrap) {
148  std::string algo = pL.get<std::string>("aggregation: drop scheme");
149 
150  TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo != "classical", Exceptions::RuntimeError, "Dropping function must not be provided for \"" << algo << "\" algorithm");
151  TEUCHOS_TEST_FOR_EXCEPTION(algo != "classical" && algo != "distance laplacian", Exceptions::RuntimeError, "\"algorithm\" must be one of (classical|distance laplacian)");
152 
153  SC threshold = as<SC>(pL.get<double>("aggregation: drop tol"));
154  GetOStream(Runtime0) << "algorithm = \"" << algo << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
155  Set(currentLevel, "Filtering", (threshold != STS::zero()));
156 
157  const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<double>("aggregation: Dirichlet threshold")));
158 
159  GO numDropped = 0, numTotal = 0;
160  std::string graphType = "unamalgamated"; //for description purposes only
161  if (algo == "classical") {
162  if (predrop_ == null) {
163  // ap: this is a hack: had to declare predrop_ as mutable
164  predrop_ = rcp(new PreDropFunctionConstVal(threshold));
165  }
166 
167  if (predrop_ != null) {
168  RCP<PreDropFunctionConstVal> predropConstVal = rcp_dynamic_cast<PreDropFunctionConstVal>(predrop_);
169  TEUCHOS_TEST_FOR_EXCEPTION(predropConstVal == Teuchos::null, Exceptions::BadCast,
170  "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
171  // If a user provided a predrop function, it overwrites the XML threshold parameter
172  SC newt = predropConstVal->GetThreshold();
173  if (newt != threshold) {
174  GetOStream(Warnings0) << "switching threshold parameter from " << threshold << " (list) to " << newt << " (user function" << std::endl;
175  threshold = newt;
176  }
177  }
178  // At this points we either have
179  // (predrop_ != null)
180  // Therefore, it is sufficient to check only threshold
181  if (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && A->hasCrsGraph()) {
182  // Case 1: scalar problem, no dropping => just use matrix graph
183  RCP<GraphBase> graph = rcp(new Graph(A->getCrsGraph(), "graph of A"));
184  // Detect and record rows that correspond to Dirichlet boundary conditions
185  ArrayRCP<const bool > boundaryNodes;
186  boundaryNodes = MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold);
187  graph->SetBoundaryNodeMap(boundaryNodes);
188  numTotal = A->getNodeNumEntries();
189 
190  if (GetVerbLevel() & Statistics1) {
191  GO numLocalBoundaryNodes = 0;
192  GO numGlobalBoundaryNodes = 0;
193  for (LO i = 0; i < boundaryNodes.size(); ++i)
194  if (boundaryNodes[i])
195  numLocalBoundaryNodes++;
196  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
197  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
198  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
199  }
200 
201  Set(currentLevel, "DofsPerNode", 1);
202  Set(currentLevel, "Graph", graph);
203 
204  } else if ( (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) ||
205  (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !A->hasCrsGraph())) {
206  // Case 2: scalar problem with dropping => record the column indices of undropped entries, but still use original
207  // graph's map information, e.g., whether index is local
208  // OR a matrix without a CrsGraph
209 
210  // allocate space for the local graph
211  ArrayRCP<LO> rows (A->getNodeNumRows()+1);
212  ArrayRCP<LO> columns(A->getNodeNumEntries());
213 
214  RCP<Vector> ghostedDiag = MueLu::Utilities<SC,LO,GO,NO>::GetMatrixOverlappedDiagonal(*A);
215  const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
216  const ArrayRCP<bool> boundaryNodes(A->getNodeNumRows(), false);
217 
218  LO realnnz = 0;
219 
220  rows[0] = 0;
221  for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
222  size_t nnz = A->getNumEntriesInLocalRow(row);
223  ArrayView<const LO> indices;
224  ArrayView<const SC> vals;
225  A->getLocalRowView(row, indices, vals);
226 
227  //FIXME the current predrop function uses the following
228  //FIXME if(std::abs(vals[k]) > std::abs(threshold_) || grow == gcid )
229  //FIXME but the threshold doesn't take into account the rows' diagonal entries
230  //FIXME For now, hardwiring the dropping in here
231 
232  LO rownnz = 0;
233  for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
234  LO col = indices[colID];
235 
236  // we avoid a square root by using squared values
237  typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]); // eps^2*|a_ii|*|a_jj|
238  typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]); // |a_ij|^2
239 
240  if (aij > aiiajj || row == col) {
241  columns[realnnz++] = col;
242  rownnz++;
243  } else
244  numDropped++;
245  }
246  if (rownnz == 1) {
247  // If the only element remaining after filtering is diagonal, mark node as boundary
248  // FIXME: this should really be replaced by the following
249  // if (indices.size() == 1 && indices[0] == row)
250  // boundaryNodes[row] = true;
251  // We do not do it this way now because there is no framework for distinguishing isolated
252  // and boundary nodes in the aggregation algorithms
253  boundaryNodes[row] = true;
254  }
255  rows[row+1] = realnnz;
256  }
257  columns.resize(realnnz);
258  numTotal = A->getNodeNumEntries();
259 
260  RCP<GraphBase> graph = rcp(new LWGraph(rows, columns, A->getRowMap(), A->getColMap(), "thresholded graph of A"));
261  graph->SetBoundaryNodeMap(boundaryNodes);
262  if (GetVerbLevel() & Statistics1) {
263  GO numLocalBoundaryNodes = 0;
264  GO numGlobalBoundaryNodes = 0;
265  for (LO i = 0; i < boundaryNodes.size(); ++i)
266  if (boundaryNodes[i])
267  numLocalBoundaryNodes++;
268  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
269  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
270  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
271  }
272  Set(currentLevel, "Graph", graph);
273  Set(currentLevel, "DofsPerNode", 1);
274 
275  } else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
276  // Case 3: Multiple DOF/node problem without dropping
277  const RCP<const Map> rowMap = A->getRowMap();
278  const RCP<const Map> colMap = A->getColMap();
279 
280  graphType = "amalgamated";
281 
282  // build node row map (uniqueMap) and node column map (nonUniqueMap)
283  // the arrays rowTranslation and colTranslation contain the local node id
284  // given a local dof id. The data is calculated by the AmalgamationFactory and
285  // stored in the variable container "UnAmalgamationInfo"
286  RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
287  RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
288  Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
289  Array<LO> colTranslation = *(amalInfo->getColTranslation());
290 
291  // get number of local nodes
292  LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
293 
294  // Allocate space for the local graph
295  ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
296  ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
297 
298  const ArrayRCP<bool> amalgBoundaryNodes(numRows, false);
299 
300  // Detect and record rows that correspond to Dirichlet boundary conditions
301  // TODO If we use ArrayRCP<LO>, then we can record boundary nodes as usual. Size
302  // TODO the array one bigger than the number of local rows, and the last entry can
303  // TODO hold the actual number of boundary nodes. Clever, huh?
304  ArrayRCP<const bool > pointBoundaryNodes;
305  pointBoundaryNodes = MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold);
306 
307  // extract striding information
308  LO blkSize = A->GetFixedBlockSize(); //< the full block size (number of dofs per node in strided map)
309  LO blkId = -1; //< the block id within the strided map (or -1 if it is a full block map)
310  LO blkPartSize = A->GetFixedBlockSize(); //< stores the size of the block within the strided map
311  if (A->IsView("stridedMaps") == true) {
312  Teuchos::RCP<const Map> myMap = A->getRowMap("stridedMaps");
313  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
314  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
315  blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
316  blkId = strMap->getStridedBlockId();
317  if (blkId > -1)
318  blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
319  }
320 
321  // loop over all local nodes
322  LO realnnz = 0;
323  rows[0] = 0;
324  Array<LO> indicesExtra;
325  for (LO row = 0; row < numRows; row++) {
326  ArrayView<const LO> indices;
327  indicesExtra.resize(0);
328 
329  // The amalgamated row is marked as Dirichlet iff all point rows are Dirichlet
330  // Note, that pointBoundaryNodes lives on the dofmap (and not the node map).
331  // Therefore, looping over all dofs is fine here. We use blkPartSize as we work
332  // with local ids.
333  // TODO: Here we have different options of how to define a node to be a boundary (or Dirichlet)
334  // node.
335  bool isBoundary = false;
336  isBoundary = true;
337  for (LO j = 0; j < blkPartSize; j++) {
338  if (!pointBoundaryNodes[row*blkPartSize+j]) {
339  isBoundary = false;
340  break;
341  }
342  }
343 
344  // Merge rows of A
345  // The array indicesExtra contains local column node ids for the current local node "row"
346  if (!isBoundary)
347  MergeRows(*A, row, indicesExtra, colTranslation);
348  else
349  indicesExtra.push_back(row);
350  indices = indicesExtra;
351  numTotal += indices.size();
352 
353  // add the local column node ids to the full columns array which
354  // contains the local column node ids for all local node rows
355  LO nnz = indices.size(), rownnz = 0;
356  for (LO colID = 0; colID < nnz; colID++) {
357  LO col = indices[colID];
358  columns[realnnz++] = col;
359  rownnz++;
360  }
361 
362  if (rownnz == 1) {
363  // If the only element remaining after filtering is diagonal, mark node as boundary
364  // FIXME: this should really be replaced by the following
365  // if (indices.size() == 1 && indices[0] == row)
366  // boundaryNodes[row] = true;
367  // We do not do it this way now because there is no framework for distinguishing isolated
368  // and boundary nodes in the aggregation algorithms
369  amalgBoundaryNodes[row] = true;
370  }
371  rows[row+1] = realnnz;
372  } //for (LO row = 0; row < numRows; row++)
373  columns.resize(realnnz);
374 
375  RCP<GraphBase> graph = rcp(new LWGraph(rows, columns, uniqueMap, nonUniqueMap, "amalgamated graph of A"));
376  graph->SetBoundaryNodeMap(amalgBoundaryNodes);
377 
378  if (GetVerbLevel() & Statistics1) {
379  GO numLocalBoundaryNodes = 0;
380  GO numGlobalBoundaryNodes = 0;
381 
382  for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
383  if (amalgBoundaryNodes[i])
384  numLocalBoundaryNodes++;
385 
386  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
387  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
388  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes
389  << " agglomerated Dirichlet nodes" << std::endl;
390  }
391 
392  Set(currentLevel, "Graph", graph);
393  Set(currentLevel, "DofsPerNode", blkSize); // full block size
394 
395  } else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
396  // Case 4: Multiple DOF/node problem with dropping
397  const RCP<const Map> rowMap = A->getRowMap();
398  const RCP<const Map> colMap = A->getColMap();
399 
400  graphType = "amalgamated";
401 
402  // build node row map (uniqueMap) and node column map (nonUniqueMap)
403  // the arrays rowTranslation and colTranslation contain the local node id
404  // given a local dof id. The data is calculated by the AmalgamationFactory and
405  // stored in the variable container "UnAmalgamationInfo"
406  RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
407  RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
408  Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
409  Array<LO> colTranslation = *(amalInfo->getColTranslation());
410 
411  // get number of local nodes
412  LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
413 
414  // Allocate space for the local graph
415  ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
416  ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
417 
418  const ArrayRCP<bool> amalgBoundaryNodes(numRows, false);
419 
420  // Detect and record rows that correspond to Dirichlet boundary conditions
421  // TODO If we use ArrayRCP<LO>, then we can record boundary nodes as usual. Size
422  // TODO the array one bigger than the number of local rows, and the last entry can
423  // TODO hold the actual number of boundary nodes. Clever, huh?
424  ArrayRCP<const bool > pointBoundaryNodes;
425  pointBoundaryNodes = MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold);
426 
427  // extract striding information
428  LO blkSize = A->GetFixedBlockSize(); //< the full block size (number of dofs per node in strided map)
429  LO blkId = -1; //< the block id within the strided map (or -1 if it is a full block map)
430  LO blkPartSize = A->GetFixedBlockSize(); //< stores the size of the block within the strided map
431  if (A->IsView("stridedMaps") == true) {
432  Teuchos::RCP<const Map> myMap = A->getRowMap("stridedMaps");
433  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
434  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
435  blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
436  blkId = strMap->getStridedBlockId();
437  if (blkId > -1)
438  blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
439  }
440 
441  // extract diagonal data for dropping strategy
442  RCP<Vector> ghostedDiag = MueLu::Utilities<SC,LO,GO,NO>::GetMatrixOverlappedDiagonal(*A);
443  const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
444 
445  // loop over all local nodes
446  LO realnnz = 0;
447  rows[0] = 0;
448  Array<LO> indicesExtra;
449  for (LO row = 0; row < numRows; row++) {
450  ArrayView<const LO> indices;
451  indicesExtra.resize(0);
452 
453  // The amalgamated row is marked as Dirichlet iff all point rows are Dirichlet
454  // Note, that pointBoundaryNodes lives on the dofmap (and not the node map).
455  // Therefore, looping over all dofs is fine here. We use blkPartSize as we work
456  // with local ids.
457  // TODO: Here we have different options of how to define a node to be a boundary (or Dirichlet)
458  // node.
459  bool isBoundary = false;
460  isBoundary = true;
461  for (LO j = 0; j < blkPartSize; j++) {
462  if (!pointBoundaryNodes[row*blkPartSize+j]) {
463  isBoundary = false;
464  break;
465  }
466  }
467 
468  // Merge rows of A
469  // The array indicesExtra contains local column node ids for the current local node "row"
470  if (!isBoundary)
471  MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
472  else
473  indicesExtra.push_back(row);
474  indices = indicesExtra;
475  numTotal += indices.size();
476 
477  // add the local column node ids to the full columns array which
478  // contains the local column node ids for all local node rows
479  LO nnz = indices.size(), rownnz = 0;
480  for (LO colID = 0; colID < nnz; colID++) {
481  LO col = indices[colID];
482  columns[realnnz++] = col;
483  rownnz++;
484  }
485 
486  if (rownnz == 1) {
487  // If the only element remaining after filtering is diagonal, mark node as boundary
488  // FIXME: this should really be replaced by the following
489  // if (indices.size() == 1 && indices[0] == row)
490  // boundaryNodes[row] = true;
491  // We do not do it this way now because there is no framework for distinguishing isolated
492  // and boundary nodes in the aggregation algorithms
493  amalgBoundaryNodes[row] = true;
494  }
495  rows[row+1] = realnnz;
496  } //for (LO row = 0; row < numRows; row++)
497  columns.resize(realnnz);
498 
499  RCP<GraphBase> graph = rcp(new LWGraph(rows, columns, uniqueMap, nonUniqueMap, "amalgamated graph of A"));
500  graph->SetBoundaryNodeMap(amalgBoundaryNodes);
501 
502  if (GetVerbLevel() & Statistics1) {
503  GO numLocalBoundaryNodes = 0;
504  GO numGlobalBoundaryNodes = 0;
505 
506  for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
507  if (amalgBoundaryNodes[i])
508  numLocalBoundaryNodes++;
509 
510  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
511  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
512  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes
513  << " agglomerated Dirichlet nodes" << std::endl;
514  }
515 
516  Set(currentLevel, "Graph", graph);
517  Set(currentLevel, "DofsPerNode", blkSize); // full block size
518  }
519 
520  } else if (algo == "distance laplacian") {
521  LO blkSize = A->GetFixedBlockSize();
522  GO indexBase = A->getRowMap()->getIndexBase();
523 
524  // [*0*] : FIXME
525  // ap: somehow, if I move this line to [*1*], Belos throws an error
526  // I'm not sure what's going on. Do we always have to Get data, if we did
527  // DeclareInput for it?
528  RCP<RealValuedMultiVector> Coords = Get< RCP<RealValuedMultiVector > >(currentLevel, "Coordinates");
529 
530  // Detect and record rows that correspond to Dirichlet boundary conditions
531  // TODO If we use ArrayRCP<LO>, then we can record boundary nodes as usual. Size
532  // TODO the array one bigger than the number of local rows, and the last entry can
533  // TODO hold the actual number of boundary nodes. Clever, huh?
534  ArrayRCP<const bool > pointBoundaryNodes;
535  pointBoundaryNodes = MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold);
536 
537  if ( (blkSize == 1) && (threshold == STS::zero()) ) {
538  // Trivial case: scalar problem, no dropping. Can return original graph
539  RCP<GraphBase> graph = rcp(new Graph(A->getCrsGraph(), "graph of A"));
540  graph->SetBoundaryNodeMap(pointBoundaryNodes);
541  graphType="unamalgamated";
542  numTotal = A->getNodeNumEntries();
543 
544  if (GetVerbLevel() & Statistics1) {
545  GO numLocalBoundaryNodes = 0;
546  GO numGlobalBoundaryNodes = 0;
547  for (LO i = 0; i < pointBoundaryNodes.size(); ++i)
548  if (pointBoundaryNodes[i])
549  numLocalBoundaryNodes++;
550  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
551  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
552  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
553  }
554 
555  Set(currentLevel, "DofsPerNode", blkSize);
556  Set(currentLevel, "Graph", graph);
557 
558  } else {
559  // ap: We make quite a few assumptions here; general case may be a lot different,
560  // but much much harder to implement. We assume that:
561  // 1) all maps are standard maps, not strided maps
562  // 2) global indices of dofs in A are related to dofs in coordinates in a simple arithmetic
563  // way: rows i*blkSize, i*blkSize+1, ..., i*blkSize + (blkSize-1) correspond to node i
564  //
565  // NOTE: Potentially, some of the code below could be simplified with UnAmalgamationInfo,
566  // but as I totally don't understand that code, here is my solution
567 
568  // [*1*]: see [*0*]
569 
570  // Check that the number of local coordinates is consistent with the #rows in A
571  TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements()/blkSize != Coords->getLocalLength(), Exceptions::Incompatible,
572  "Coordinate vector length (" << Coords->getLocalLength() << ") is incompatible with number of rows in A (" << A->getRowMap()->getNodeNumElements() << ") by modulo block size ("<< blkSize <<").");
573 
574  const RCP<const Map> colMap = A->getColMap();
575  RCP<const Map> uniqueMap, nonUniqueMap;
576  Array<LO> colTranslation;
577  if (blkSize == 1) {
578  uniqueMap = A->getRowMap();
579  nonUniqueMap = A->getColMap();
580  graphType="unamalgamated";
581 
582  } else {
583  uniqueMap = Coords->getMap();
584  TEUCHOS_TEST_FOR_EXCEPTION(uniqueMap->getIndexBase() != indexBase, Exceptions::Incompatible,
585  "Different index bases for matrix and coordinates");
586 
587  AmalgamationFactory::AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, colTranslation);
588 
589  graphType = "amalgamated";
590  }
591  LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
592 
593  RCP<RealValuedMultiVector> ghostedCoords;
594  RCP<Vector> ghostedLaplDiag;
595  Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
596  if (threshold != STS::zero()) {
597  // Get ghost coordinates
598  RCP<const Import> importer;
599  {
600  SubFactoryMonitor m1(*this, "Import construction", currentLevel);
601  if (blkSize == 1 && A->getCrsGraph()->getImporter() != Teuchos::null) {
602  GetOStream(Warnings1) << "Using existing importer from matrix graph" << std::endl;
603  importer = A->getCrsGraph()->getImporter();
604  } else {
605  GetOStream(Warnings0) << "Constructing new importer instance" << std::endl;
606  importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
607  }
608  } //subtimer
609  ghostedCoords = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(nonUniqueMap, Coords->getNumVectors());
610  {
611  SubFactoryMonitor m1(*this, "Coordinate import", currentLevel);
612  ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
613  } //subtimer
614 
615  // Construct Distance Laplacian diagonal
616  RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
617  ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
618  Array<LO> indicesExtra;
619  Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
620  if (threshold != STS::zero()) {
621  const size_t numVectors = ghostedCoords->getNumVectors();
622  coordData.reserve(numVectors);
623  for (size_t j = 0; j < numVectors; j++) {
624  Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
625  coordData.push_back(tmpData);
626  }
627  }
628  {
629  SubFactoryMonitor m1(*this, "Laplacian local diagonal", currentLevel);
630  for (LO row = 0; row < numRows; row++) {
631  ArrayView<const LO> indices;
632 
633  if (blkSize == 1) {
634  ArrayView<const SC> vals;
635  A->getLocalRowView(row, indices, vals);
636 
637  } else {
638  // Merge rows of A
639  indicesExtra.resize(0);
640  MergeRows(*A, row, indicesExtra, colTranslation);
641  indices = indicesExtra;
642  }
643 
644  LO nnz = indices.size();
645  for (LO colID = 0; colID < nnz; colID++) {
646  const LO col = indices[colID];
647 
648  if (row != col) {
649  localLaplDiagData[row] += STS::one() / MueLu::Utilities<real_type,LO,GO,NO>::Distance2(coordData, row, col);
650  }
651  }
652  }
653  } //subtimer
654  {
655  SubFactoryMonitor m1(*this, "Laplacian distributed diagonal", currentLevel);
656  ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
657  ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
658  ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
659  } //subtimer
660 
661  } else {
662  GetOStream(Runtime0) << "Skipping distance laplacian construction due to 0 threshold" << std::endl;
663  }
664 
665  // NOTE: ghostedLaplDiagData might be zero if we don't actually calculate the laplacian
666 
667  // allocate space for the local graph
668  ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
669  ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
670 
671  const ArrayRCP<bool> amalgBoundaryNodes(numRows, false);
672 
673  LO realnnz = 0;
674  rows[0] = 0;
675  Array<LO> indicesExtra;
676  {
677  SubFactoryMonitor m1(*this, "Laplacian dropping", currentLevel);
678  Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
679  if (threshold != STS::zero()) {
680  const size_t numVectors = ghostedCoords->getNumVectors();
681  coordData.reserve(numVectors);
682  for (size_t j = 0; j < numVectors; j++) {
683  Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
684  coordData.push_back(tmpData);
685  }
686  }
687  for (LO row = 0; row < numRows; row++) {
688  ArrayView<const LO> indices;
689  indicesExtra.resize(0);
690 
691  if (blkSize == 1) {
692  ArrayView<const SC> vals;
693  A->getLocalRowView(row, indices, vals);
694 
695  } else {
696  // The amalgamated row is marked as Dirichlet iff all point rows are Dirichlet
697  bool isBoundary = false;
698  isBoundary = true;
699  for (LO j = 0; j < blkSize; j++) {
700  if (!pointBoundaryNodes[row*blkSize+j]) {
701  isBoundary = false;
702  break;
703  }
704  }
705 
706  // Merge rows of A
707  if (!isBoundary)
708  MergeRows(*A, row, indicesExtra, colTranslation);
709  else
710  indicesExtra.push_back(row);
711  indices = indicesExtra;
712  }
713  numTotal += indices.size();
714 
715  LO nnz = indices.size(), rownnz = 0;
716  if (threshold != STS::zero()) {
717  for (LO colID = 0; colID < nnz; colID++) {
718  LO col = indices[colID];
719 
720  if (row == col) {
721  columns[realnnz++] = col;
722  rownnz++;
723  continue;
724  }
725 
726  SC laplVal = STS::one() / MueLu::Utilities<real_type,LO,GO,NO>::Distance2(coordData, row, col);
727  typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
728  typename STS::magnitudeType aij = STS::magnitude(laplVal*laplVal);
729 
730  if (aij > aiiajj) {
731  columns[realnnz++] = col;
732  rownnz++;
733  } else {
734  numDropped++;
735  }
736  }
737 
738  } else {
739  // Skip laplace calculation and threshold comparison for zero threshold
740  for (LO colID = 0; colID < nnz; colID++) {
741  LO col = indices[colID];
742  columns[realnnz++] = col;
743  rownnz++;
744  }
745  }
746 
747  if (rownnz == 1) {
748  // If the only element remaining after filtering is diagonal, mark node as boundary
749  // FIXME: this should really be replaced by the following
750  // if (indices.size() == 1 && indices[0] == row)
751  // boundaryNodes[row] = true;
752  // We do not do it this way now because there is no framework for distinguishing isolated
753  // and boundary nodes in the aggregation algorithms
754  amalgBoundaryNodes[row] = true;
755  }
756  rows[row+1] = realnnz;
757  } //for (LO row = 0; row < numRows; row++)
758  } //subtimer
759  columns.resize(realnnz);
760 
761  RCP<GraphBase> graph;
762  {
763  SubFactoryMonitor m1(*this, "Build amalgamated graph", currentLevel);
764  graph = rcp(new LWGraph(rows, columns, uniqueMap, nonUniqueMap, "amalgamated graph of A"));
765  graph->SetBoundaryNodeMap(amalgBoundaryNodes);
766  } //subtimer
767 
768  if (GetVerbLevel() & Statistics1) {
769  GO numLocalBoundaryNodes = 0;
770  GO numGlobalBoundaryNodes = 0;
771 
772  for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
773  if (amalgBoundaryNodes[i])
774  numLocalBoundaryNodes++;
775 
776  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
777  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
778  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " agglomerated Dirichlet nodes"
779  << " using threshold " << dirichletThreshold << std::endl;
780  }
781 
782  Set(currentLevel, "Graph", graph);
783  Set(currentLevel, "DofsPerNode", blkSize);
784  }
785  }
786 
787  if ((GetVerbLevel() & Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
788  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
789  GO numGlobalTotal, numGlobalDropped;
790  MueLu_sumAll(comm, numTotal, numGlobalTotal);
791  MueLu_sumAll(comm, numDropped, numGlobalDropped);
792  GetOStream(Statistics1) << "Number of dropped entries in " << graphType << " matrix graph: " << numGlobalDropped << "/" << numGlobalTotal;
793  if (numGlobalTotal != 0)
794  GetOStream(Statistics1) << " (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) << "%)";
795  GetOStream(Statistics1) << std::endl;
796  }
797 
798  } else {
799  //what Tobias has implemented
800 
801  SC threshold = as<SC>(pL.get<double>("aggregation: drop tol"));
802  //GetOStream(Runtime0) << "algorithm = \"" << algo << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
803  GetOStream(Runtime0) << "algorithm = \"" << "failsafe" << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
804  Set(currentLevel, "Filtering", (threshold != STS::zero()));
805 
806  RCP<const Map> rowMap = A->getRowMap();
807  RCP<const Map> colMap = A->getColMap();
808 
809  LO blockdim = 1; // block dim for fixed size blocks
810  GO indexBase = rowMap->getIndexBase(); // index base of maps
811  GO offset = 0;
812 
813  // 1) check for blocking/striding information
814  if(A->IsView("stridedMaps") &&
815  Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
816  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
817  RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
818  TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,Exceptions::BadCast,"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
819  blockdim = strMap->getFixedBlockSize();
820  offset = strMap->getOffset();
821  oldView = A->SwitchToView(oldView);
822  GetOStream(Statistics1) << "CoalesceDropFactory::Build():" << " found blockdim=" << blockdim << " from strided maps. offset=" << offset << std::endl;
823  } else GetOStream(Statistics1) << "CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
824 
825  // 2) get row map for amalgamated matrix (graph of A)
826  // with same distribution over all procs as row map of A
827  RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
828  GetOStream(Statistics1) << "CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() << "/" << nodeMap->getGlobalNumElements() << " elements" << std::endl;
829 
830  // 3) create graph of amalgamated matrix
831  RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getNodeMaxNumRowEntries()*blockdim, Xpetra::StaticProfile);
832 
833  LO numRows = A->getRowMap()->getNodeNumElements();
834  LO numNodes = nodeMap->getNodeNumElements();
835  const ArrayRCP<bool> amalgBoundaryNodes(numNodes, false);
836  const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0); // helper array counting the number of Dirichlet nodes associated with node
837  bool bIsDiagonalEntry = false; // boolean flag stating that grid==gcid
838 
839  // 4) do amalgamation. generate graph of amalgamated matrix
840  // Note, this code is much more inefficient than the leightwight implementation
841  // Most of the work has already been done in the AmalgamationFactory
842  for(LO row=0; row<numRows; row++) {
843  // get global DOF id
844  GO grid = rowMap->getGlobalElement(row);
845 
846  // reinitialize boolean helper variable
847  bIsDiagonalEntry = false;
848 
849  // translate grid to nodeid
850  GO nodeId = AmalgamationFactory::DOFGid2NodeId(grid, blockdim, offset, indexBase);
851 
852  size_t nnz = A->getNumEntriesInLocalRow(row);
853  Teuchos::ArrayView<const LO> indices;
854  Teuchos::ArrayView<const SC> vals;
855  A->getLocalRowView(row, indices, vals);
856 
857  RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(new std::vector<GO>); // global column block ids
858  LO realnnz = 0;
859  for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
860  GO gcid = colMap->getGlobalElement(indices[col]); // global column id
861 
862  if(vals[col]!=STS::zero()) {
863  GO cnodeId = AmalgamationFactory::DOFGid2NodeId(gcid, blockdim, offset, indexBase);
864  cnodeIds->push_back(cnodeId);
865  realnnz++; // increment number of nnz in matrix row
866  if (grid == gcid) bIsDiagonalEntry = true;
867  }
868  }
869 
870  if(realnnz == 1 && bIsDiagonalEntry == true) {
871  LO lNodeId = nodeMap->getLocalElement(nodeId);
872  numberDirichletRowsPerNode[lNodeId] += 1; // increment Dirichlet row counter associated with lNodeId
873  if (numberDirichletRowsPerNode[lNodeId] == blockdim) // mark full Dirichlet nodes
874  amalgBoundaryNodes[lNodeId] = true;
875  }
876 
877  Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
878 
879  if(arr_cnodeIds.size() > 0 )
880  crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
881  }
882  // fill matrix graph
883  crsGraph->fillComplete(nodeMap,nodeMap);
884 
885  // 5) create MueLu Graph object
886  RCP<GraphBase> graph = rcp(new Graph(crsGraph, "amalgamated graph of A"));
887 
888  // Detect and record rows that correspond to Dirichlet boundary conditions
889  graph->SetBoundaryNodeMap(amalgBoundaryNodes);
890 
891  if (GetVerbLevel() & Statistics1) {
892  GO numLocalBoundaryNodes = 0;
893  GO numGlobalBoundaryNodes = 0;
894  for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
895  if (amalgBoundaryNodes[i])
896  numLocalBoundaryNodes++;
897  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
898  MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
899  GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
900  }
901 
902  // 6) store results in Level
903  //graph->SetBoundaryNodeMap(gBoundaryNodeMap);
904  Set(currentLevel, "DofsPerNode", blockdim);
905  Set(currentLevel, "Graph", graph);
906 
907  } //if (doExperimentalWrap) ... else ...
908 
909  } //Build
910 
911  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
912  void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MergeRows(const Matrix& A, const LO row, Array<LO>& cols, const Array<LO>& translation) const {
913  typedef typename ArrayView<const LO>::size_type size_type;
914 
915  // extract striding information
916  LO blkSize = A.GetFixedBlockSize(); //< stores the size of the block within the strided map
917  if (A.IsView("stridedMaps") == true) {
918  Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
919  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
920  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
921  if (strMap->getStridedBlockId() > -1)
922  blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
923  }
924 
925  // count nonzero entries in all dof rows associated with node row
926  size_t nnz = 0, pos = 0;
927  for (LO j = 0; j < blkSize; j++)
928  nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
929 
930  if (nnz == 0) {
931  cols.resize(0);
932  return;
933  }
934 
935  cols.resize(nnz);
936 
937  // loop over all local dof rows associated with local node "row"
938  ArrayView<const LO> inds;
939  ArrayView<const SC> vals;
940  for (LO j = 0; j < blkSize; j++) {
941  A.getLocalRowView(row*blkSize+j, inds, vals);
942  size_type numIndices = inds.size();
943 
944  if (numIndices == 0) // skip empty dof rows
945  continue;
946 
947  // cols: stores all local node ids for current local node id "row"
948  cols[pos++] = translation[inds[0]];
949  for (size_type k = 1; k < numIndices; k++) {
950  LO nodeID = translation[inds[k]];
951  // Here we try to speed up the process by reducing the size of an array
952  // to sort. This works if the column nonzeros belonging to the same
953  // node are stored consequently.
954  if (nodeID != cols[pos-1])
955  cols[pos++] = nodeID;
956  }
957  }
958  cols.resize(pos);
959  nnz = pos;
960 
961  // Sort and remove duplicates
962  std::sort(cols.begin(), cols.end());
963  pos = 0;
964  for (size_t j = 1; j < nnz; j++)
965  if (cols[j] != cols[pos])
966  cols[++pos] = cols[j];
967  cols.resize(pos+1);
968  }
969 
970  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
971  void CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MergeRowsWithDropping(const Matrix& A, const LO row, const ArrayRCP<const SC>& ghostedDiagVals, SC threshold, Array<LO>& cols, const Array<LO>& translation) const {
972  typedef typename ArrayView<const LO>::size_type size_type;
973  typedef Teuchos::ScalarTraits<SC> STS;
974 
975  // extract striding information
976  LO blkSize = A.GetFixedBlockSize(); //< stores the size of the block within the strided map
977  if (A.IsView("stridedMaps") == true) {
978  Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
979  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
980  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
981  if (strMap->getStridedBlockId() > -1)
982  blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
983  }
984 
985  // count nonzero entries in all dof rows associated with node row
986  size_t nnz = 0, pos = 0;
987  for (LO j = 0; j < blkSize; j++)
988  nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
989 
990  if (nnz == 0) {
991  cols.resize(0);
992  return;
993  }
994 
995  cols.resize(nnz);
996 
997  // loop over all local dof rows associated with local node "row"
998  ArrayView<const LO> inds;
999  ArrayView<const SC> vals;
1000  for (LO j = 0; j < blkSize; j++) {
1001  A.getLocalRowView(row*blkSize+j, inds, vals);
1002  size_type numIndices = inds.size();
1003 
1004  if (numIndices == 0) // skip empty dof rows
1005  continue;
1006 
1007  // cols: stores all local node ids for current local node id "row"
1008  LO prevNodeID = -1;
1009  for (size_type k = 0; k < numIndices; k++) {
1010  LO dofID = inds[k];
1011  LO nodeID = translation[inds[k]];
1012 
1013  // we avoid a square root by using squared values
1014  typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]); // eps^2 * |a_ii| * |a_jj|
1015  typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
1016 
1017  // check dropping criterion
1018  if (aij > aiiajj || (row*blkSize+j == dofID)) {
1019  // accept entry in graph
1020 
1021  // Here we try to speed up the process by reducing the size of an array
1022  // to sort. This works if the column nonzeros belonging to the same
1023  // node are stored consequently.
1024  if (nodeID != prevNodeID) {
1025  cols[pos++] = nodeID;
1026  prevNodeID = nodeID;
1027  }
1028  }
1029  }
1030  }
1031  cols.resize(pos);
1032  nnz = pos;
1033 
1034  // Sort and remove duplicates
1035  std::sort(cols.begin(), cols.end());
1036  pos = 0;
1037  for (size_t j = 1; j < nnz; j++)
1038  if (cols[j] != cols[pos])
1039  cols[++pos] = cols[j];
1040  cols.resize(pos+1);
1041 
1042  return;
1043  }
1044 } //namespace MueLu
1045 
1046 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP
MueLu::Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:57
MueLu::CoalesceDropFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_CoalesceDropFactory_def.hpp:79
MueLu::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: MueLu_Exceptions.hpp:57
Teuchos::Array::reserve
void reserve(size_type n)
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:61
Xpetra_Vector.hpp
MueLu_CoalesceDropFactory_decl.hpp
Teuchos::ParameterEntry::setValidator
void setValidator(RCP< const ParameterEntryValidator > const &validator)
Xpetra::MultiVector::description
virtual std::string description() const =0
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu_Exceptions.hpp
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::Graph
MueLu representation of a compressed row storage graph.
Definition: MueLu_Graph_decl.hpp:70
MueLu::Utilities::DetectDirichletRows
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
Definition: MueLu_Utilities_decl.hpp:219
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Xpetra_CrsGraphFactory.hpp
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
Teuchos::ArrayView::size
size_type size() const
LO
LocalOrdinal LO
MueLu::Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:62
MueLu_Level.hpp
Xpetra_Matrix.hpp
Teuchos::ArrayRCP::size
size_type size() const
MueLu::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: MueLu_Exceptions.hpp:82
MueLu::Utilities::Distance2
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Definition: MueLu_Utilities_decl.hpp:218
Xpetra_CrsGraph.hpp
MueLu::AmalgamationFactory::DOFGid2NodeId
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
Definition: MueLu_AmalgamationFactory_def.hpp:206
MueLu_MasterList.hpp
MueLu::GraphBase::SetBoundaryNodeMap
virtual void SetBoundaryNodeMap(const ArrayRCP< const bool > &boundaryArray)=0
MueLu_GraphBase.hpp
MueLu::PreDropFunctionConstVal::GetThreshold
Scalar GetThreshold() const
Return threshold value.
Definition: MueLu_PreDropFunctionConstVal_def.hpp:70
MueLu::Warnings1
Additional warnings.
Definition: MueLu_VerbosityLevel.hpp:59
Teuchos::ArrayRCP::resize
void resize(const size_type n, const T &val=T())
Xpetra_MultiVectorFactory.hpp
Xpetra_MultiVector.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:261
Xpetra_StridedMap.hpp
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
MueLu::AmalgamationFactory::AmalgamateMap
static void AmalgamateMap(const Map &sourceMap, const Matrix &A, RCP< const Map > &amalgamatedMap, Array< LO > &translation)
Method to create merged map for systems of PDEs.
Definition: MueLu_AmalgamationFactory_def.hpp:159
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
MueLu::CoalesceDropFactory::Build
void Build(Level &currentLevel) const
Build an object with this factory.
Definition: MueLu_CoalesceDropFactory_def.hpp:118
Xpetra::viewLabel_t
std::string viewLabel_t
Teuchos::Array::end
iterator end()
MueLu::CoalesceDropFactory::DeclareInput
void DeclareInput(Level &currentLevel) const
Input.
Definition: MueLu_CoalesceDropFactory_def.hpp:105
MueLu::Parameters0
Print class parameters.
Definition: MueLu_VerbosityLevel.hpp:67
Xpetra_MapFactory.hpp
Teuchos::Array::push_back
void push_back(const value_type &x)
MueLu::CoalesceDropFactory::MergeRows
void MergeRows(const Matrix &A, const LO row, Array< LO > &cols, const Array< LO > &translation) const
Method to merge rows of matrix for systems of PDEs.
Definition: MueLu_CoalesceDropFactory_def.hpp:912
MueLu::Utilities::GetMatrixOverlappedDiagonal
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Definition: MueLu_Utilities_decl.hpp:210
SC
Scalar SC
MueLu_Monitor.hpp
MueLu::LWGraph
Lightweight MueLu representation of a compressed row storage graph.
Definition: MueLu_LWGraph_decl.hpp:70
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Xpetra_VectorFactory.hpp
Xpetra_Map.hpp
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
Teuchos::ParameterList::getEntry
ParameterEntry & getEntry(const std::string &name)
Teuchos::Array::begin
iterator begin()
MueLu::CoalesceDropFactory::MergeRowsWithDropping
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const
Definition: MueLu_CoalesceDropFactory_def.hpp:971
Xpetra_ImportFactory.hpp