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