MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_AmalgamationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_AMALGAMATIONFACTORY_DEF_HPP
11 #define MUELU_AMALGAMATIONFACTORY_DEF_HPP
12 
13 #include <Xpetra_Matrix.hpp>
14 
16 
17 #include "MueLu_Level.hpp"
18 #include "MueLu_AmalgamationInfo.hpp"
19 #include "MueLu_Monitor.hpp"
20 
21 namespace MueLu {
22 
23 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
25 
26 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
28 
29 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
31  RCP<ParameterList> validParamList = rcp(new ParameterList());
32  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
33  return validParamList;
34 }
35 
36 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
38  Input(currentLevel, "A"); // sub-block from blocked A
39 }
40 
41 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43  FactoryMonitor m(*this, "Build", currentLevel);
44 
45  RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel, "A");
46 
47  /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
48  fullblocksize is the number of storage blocks that must kept together during the amalgamation process.
49 
50  Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
51 
52  numPDEs = fullblocksize * storageblocksize.
53 
54  If numPDEs==1
55  Matrix is point storage (classical CRS storage). storageblocksize=1 and fullblocksize=1
56  No other values makes sense.
57 
58  If numPDEs>1
59  If matrix uses point storage, then storageblocksize=1 and fullblockssize=numPDEs.
60  If matrix uses block storage, with block size of n, then storageblocksize=n, and fullblocksize=numPDEs/n.
61  Thus far, only storageblocksize=numPDEs and fullblocksize=1 has been tested.
62  */
63 
64  LO fullblocksize = 1; // block dim for fixed size blocks
65  GO offset = 0; // global offset of dof gids
66  LO blockid = -1; // block id in strided map
67  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
68  LO stridedblocksize = fullblocksize; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
69  LO storageblocksize = A->GetStorageBlockSize();
70  // GO indexBase = A->getRowMap()->getIndexBase(); // index base for maps (unused)
71 
72  // 1) check for blocking/striding information
73 
74  if (A->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
75  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // NOTE: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
76  RCP<const StridedMap> stridedRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
77  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap == Teuchos::null, Exceptions::BadCast, "MueLu::CoalesceFactory::Build: cast to strided row map failed.");
78  fullblocksize = stridedRowMap->getFixedBlockSize();
79  offset = DOFGidOffset(stridedRowMap);
80  blockid = stridedRowMap->getStridedBlockId();
81 
82  if (blockid > -1) {
83  std::vector<size_t> stridingInfo = stridedRowMap->getStridingData();
84  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
85  nStridedOffset += stridingInfo[j];
86  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
87 
88  } else {
89  stridedblocksize = fullblocksize;
90  }
91  // Correct for the storageblocksize
92  // NOTE: Before this point fullblocksize is actually numPDEs
93  TEUCHOS_TEST_FOR_EXCEPTION(fullblocksize % storageblocksize != 0, Exceptions::RuntimeError, "AmalgamationFactory: fullblocksize needs to be a multiple of A->GetStorageBlockSize()");
94  fullblocksize /= storageblocksize;
95  stridedblocksize /= storageblocksize;
96 
97  oldView = A->SwitchToView(oldView);
98  GetOStream(Runtime1) << "AmalagamationFactory::Build():"
99  << " found fullblocksize=" << fullblocksize << " and stridedblocksize=" << stridedblocksize << " from strided maps. offset=" << offset << std::endl;
100 
101  } else {
102  GetOStream(Warnings0) << "AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
103  }
104 
105  // build node row map (uniqueMap) and node column map (nonUniqueMap)
106  // the arrays rowTranslation and colTranslation contain the local node id
107  // given a local dof id. They are only necessary for the CoalesceDropFactory if
108  // fullblocksize > 1
109  RCP<const Map> uniqueMap, nonUniqueMap;
110  RCP<AmalgamationInfo> amalgamationData;
111  RCP<Array<LO> > rowTranslation = Teuchos::null;
112  RCP<Array<LO> > colTranslation = Teuchos::null;
113 
114  if (fullblocksize > 1) {
115  // mfh 14 Apr 2015: These need to have different names than
116  // rowTranslation and colTranslation, in order to avoid
117  // shadowing warnings (-Wshadow with GCC). Alternately, it
118  // looks like you could just assign to the existing variables in
119  // this scope, rather than creating new ones.
120  RCP<Array<LO> > theRowTranslation = rcp(new Array<LO>);
121  RCP<Array<LO> > theColTranslation = rcp(new Array<LO>);
122  AmalgamateMap(*(A->getRowMap()), *A, uniqueMap, *theRowTranslation);
123  AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, *theColTranslation);
124 
125  amalgamationData = rcp(new AmalgamationInfo(theRowTranslation,
126  theColTranslation,
127  uniqueMap,
128  nonUniqueMap,
129  A->getColMap(),
130  fullblocksize,
131  offset,
132  blockid,
133  nStridedOffset,
134  stridedblocksize));
135  } else {
136  amalgamationData = rcp(new AmalgamationInfo(rowTranslation, // Teuchos::null
137  colTranslation, // Teuchos::null
138  A->getRowMap(), // unique map of graph
139  A->getColMap(), // non-unique map of graph
140  A->getColMap(), // column map of A
141  fullblocksize,
142  offset,
143  blockid,
144  nStridedOffset,
145  stridedblocksize));
146  }
147 
148  // store (un)amalgamation information on current level
149  Set(currentLevel, "UnAmalgamationInfo", amalgamationData);
150 }
151 
152 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
153 void AmalgamationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AmalgamateMap(const Map& sourceMap, const Matrix& A, RCP<const Map>& amalgamatedMap, Array<LO>& translation) {
154  typedef typename ArrayView<const GO>::size_type size_type;
155  typedef std::unordered_map<GO, size_type> container;
156 
157  GO indexBase = sourceMap.getIndexBase();
158  ArrayView<const GO> elementAList = sourceMap.getLocalElementList();
159  size_type numElements = elementAList.size();
160  container filter;
161 
162  GO offset = 0;
163  LO blkSize = A.GetFixedBlockSize() / A.GetStorageBlockSize();
164  if (A.IsView("stridedMaps") == true) {
165  Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
166  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
167  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
168  offset = DOFGidOffset(strMap);
169  blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
170  }
171 
172  Array<GO> elementList(numElements);
173  translation.resize(numElements);
174 
175  size_type numRows = 0;
176  for (size_type id = 0; id < numElements; id++) {
177  GO dofID = elementAList[id];
178  GO nodeID = AmalgamationFactory::DOFGid2NodeId(dofID, blkSize, offset, indexBase);
179 
180  typename container::iterator it = filter.find(nodeID);
181  if (it == filter.end()) {
182  filter[nodeID] = numRows;
183 
184  translation[id] = numRows;
185  elementList[numRows] = nodeID;
186 
187  numRows++;
188 
189  } else {
190  translation[id] = it->second;
191  }
192  }
193  elementList.resize(numRows);
194 
195  amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
196 }
197 
198 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
200  // NOTE: There could be further optimizations here where we detect contiguous maps and then
201  // create a contiguous amalgamated maps, which bypasses the expense of the getMyGlobalIndicesDevice()
202  // call (which is free for non-contiguous maps, but costs something if the map is contiguous).
203  using range_policy = Kokkos::RangePolicy<typename Node::execution_space>;
204  using array_type = typename Map::global_indices_array_device_type;
205 
206  array_type elementAList = sourceMap->getMyGlobalIndicesDevice();
207  GO indexBase = sourceMap->getIndexBase();
208  LO blkSize = sourceMap->getFixedBlockSize();
209  auto numElements = elementAList.size() / blkSize;
210  typename array_type::non_const_type elementList_nc("elementList", numElements);
211 
212  // Amalgamate the map
213  Kokkos::parallel_for(
214  "Amalgamate Element List", range_policy(0, numElements), KOKKOS_LAMBDA(LO i) {
215  elementList_nc[i] = (elementAList[i * blkSize] - indexBase) / blkSize + indexBase;
216  });
217  array_type elementList = elementList_nc;
218 
219  amalgamatedMap = MapFactory::Build(sourceMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
220  elementList, indexBase, sourceMap->getComm());
221 }
222 
223 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
225  const GlobalOrdinal offset, const GlobalOrdinal indexBase) {
226  GlobalOrdinal globalblockid = ((GlobalOrdinal)gid - offset - indexBase) / blockSize + indexBase;
227  return globalblockid;
228 }
229 
230 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
232  GlobalOrdinal offset = stridedRowMap->getMinAllGlobalIndex();
233 
234  size_t nStridedOffset = 0;
235  for (int j = 0; j < stridedRowMap->getStridedBlockId(); j++) {
236  nStridedOffset += stridedRowMap->getStridingData()[j];
237  }
238  const GO goStridedOffset = Teuchos::as<const GO>(nStridedOffset);
239 
240  offset -= goStridedOffset + stridedRowMap->getIndexBase();
241 
242  return offset;
243 }
244 
245 } // namespace MueLu
246 
247 #endif /* MUELU_SUBBLOCKUNAMALGAMATIONFACTORY_DEF_HPP */
Important warning messages (one line)
void DeclareInput(Level &currentLevel) const override
Input.
Exception indicating invalid cast attempted.
MueLu::DefaultLocalOrdinal LocalOrdinal
GlobalOrdinal GO
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
LocalOrdinal LO
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static const GlobalOrdinal DOFGidOffset(RCP< const StridedMap > stridedMap)
Method to calculate the global (row) id offset from scratch.
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
AmalgamationFactory()
Constructor.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
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 resize(size_type new_size, const value_type &x=value_type())
void Build(Level &currentLevel) const override
Build an object with this factory.
std::string viewLabel_t
virtual ~AmalgamationFactory()
Destructor.
RCP< const ParameterList > GetValidParameterList() const override
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
minimal container class for storing amalgamation information