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 
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>
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 
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() & Statistics1) {
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(Statistics1) << "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() & Statistics1) {
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(Statistics1) << "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() & Statistics1) {
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(Statistics1) << "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() & Statistics1) {
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(Statistics1) << "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() & Statistics1) {
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(Statistics1) << "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() & Statistics1) {
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(Statistics1) << "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() & Statistics1) && !(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(Statistics1) << "Number of dropped entries in " << graphType << " matrix graph: " << numGlobalDropped << "/" << numGlobalTotal;
756  if (numGlobalTotal != 0)
757  GetOStream(Statistics1) << " (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) << "%)";
758  GetOStream(Statistics1) << 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(Statistics1) << "CoalesceDropFactory::Build():" << " found blockdim=" << blockdim << " from strided maps. offset=" << offset << std::endl;
787  } else GetOStream(Statistics1) << "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(Statistics1) << "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() & Statistics1) {
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(Statistics1) << "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)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
virtual std::string description() const=0
GlobalOrdinal GO
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.
Print more statistics.
LocalOrdinal LO
One-liner description of what is happening.
Namespace for MueLu classes and methods.
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
void DeclareInput(Level &currentLevel) const
Input.
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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#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 MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const
void Build(Level &currentLevel) const
Build an object with this factory.
std::string viewLabel_t
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)
Scalar SC
Lightweight MueLu representation of a compressed row storage graph.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
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.