MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_AmalgamationFactory_kokkos_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_AMALGAMATIONFACTORY_KOKKOS_DEF_HPP
47 #define MUELU_AMALGAMATIONFACTORY_KOKKOS_DEF_HPP
48 
49 #include <Xpetra_Matrix.hpp>
50 
51 
52 #include "MueLu_Level.hpp"
53 #include "MueLu_Monitor.hpp"
54 
56 #include "MueLu_AmalgamationInfo_kokkos.hpp"
57 
58 namespace MueLu {
59 
60  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63  RCP<ParameterList> validParamList = rcp(new ParameterList());
64  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
65  return validParamList;
66  }
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  DeclareInput(Level &currentLevel) const {
71  Input(currentLevel, "A"); // sub-block from blocked A
72  }
73 
74  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
76  Build(Level &currentLevel) const
77  {
78  FactoryMonitor m(*this, "Build", currentLevel);
79 
80  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
81 
82  /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
83  fullblocksize is the number of storage blocks that must kept together during the amalgamation process.
84 
85  Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
86 
87  numPDEs = fullblocksize * storageblocksize.
88 
89  If numPDEs==1
90  Matrix is point storage (classical CRS storage). storageblocksize=1 and fullblocksize=1
91  No other values makes sense.
92 
93  If numPDEs>1
94  If matrix uses point storage, then storageblocksize=1 and fullblockssize=numPDEs.
95  If matrix uses block storage, with block size of n, then storageblocksize=n, and fullblocksize=numPDEs/n.
96  Thus far, only storageblocksize=numPDEs and fullblocksize=1 has been tested.
97  */
98 
99  LO fullblocksize = 1; // block dim for fixed size blocks
100  GO offset = 0; // global offset of dof gids
101  LO blockid = -1; // block id in strided map
102  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
103  LO stridedblocksize = fullblocksize; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
104  LO storageblocksize = A->GetStorageBlockSize();
105  // GO indexBase = A->getRowMap()->getIndexBase(); // index base for maps (unused)
106 
107  // 1) check for blocking/striding information
108 
109  if (A->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
110  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // NOTE: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
111  RCP<const StridedMap> stridedRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
112  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap == Teuchos::null,Exceptions::BadCast,"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
113  fullblocksize = stridedRowMap->getFixedBlockSize();
114  offset = stridedRowMap->getOffset();
115  blockid = stridedRowMap->getStridedBlockId();
116 
117  if (blockid > -1) {
118  std::vector<size_t> stridingInfo = stridedRowMap->getStridingData();
119  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
120  nStridedOffset += stridingInfo[j];
121  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
122 
123  } else {
124  stridedblocksize = fullblocksize;
125  }
126  // Correct for the storageblocksize
127  // NOTE: Before this point fullblocksize is actually numPDEs
128  TEUCHOS_TEST_FOR_EXCEPTION(fullblocksize % storageblocksize != 0,Exceptions::RuntimeError,"AmalgamationFactory::Build(): fullblocksize needs to be a multiple of A->GetStorageBlockSize()");
129  fullblocksize /= storageblocksize;
130  stridedblocksize /= storageblocksize;
131 
132  oldView = A->SwitchToView(oldView);
133  GetOStream(Runtime1) << "AmalagamationFactory::Build():" << " found fullblocksize=" << fullblocksize << " and stridedblocksize=" << stridedblocksize << " from strided maps. offset=" << offset << std::endl;
134 
135  } else {
136  GetOStream(Warnings0) << "AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
137  }
138 
139  // build node row map (uniqueMap) and node column map (nonUniqueMap)
140  // the arrays rowTranslation and colTranslation contain the local node id
141  // given a local dof id. They are only necessary for the CoalesceDropFactory if
142  // fullblocksize > 1
143  RCP<const Map> uniqueMap, nonUniqueMap;
144  RCP<AmalgamationInfo_kokkos> amalgamationData;
145  RCP<Array<LO> > rowTranslation = Teuchos::null;
146  RCP<Array<LO> > colTranslation = Teuchos::null;
147 
148  if (fullblocksize > 1) {
149  // mfh 14 Apr 2015: These need to have different names than
150  // rowTranslation and colTranslation, in order to avoid
151  // shadowing warnings (-Wshadow with GCC). Alternately, it
152  // looks like you could just assign to the existing variables in
153  // this scope, rather than creating new ones.
154  RCP<Array<LO> > theRowTranslation = rcp(new Array<LO>);
155  RCP<Array<LO> > theColTranslation = rcp(new Array<LO>);
156  AmalgamateMap(*(A->getRowMap()), *A, uniqueMap, *theRowTranslation);
157  AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, *theColTranslation);
158 
159  amalgamationData = rcp(new AmalgamationInfo_kokkos(theRowTranslation,
160  theColTranslation,
161  uniqueMap,
162  nonUniqueMap,
163  A->getColMap(),
164  fullblocksize,
165  offset,
166  blockid,
167  nStridedOffset,
168  stridedblocksize) );
169  } else {
170  amalgamationData = rcp(new AmalgamationInfo_kokkos(rowTranslation, // Teuchos::null
171  colTranslation, // Teuchos::null
172  A->getRowMap(), // unique map of graph
173  A->getColMap(), // non-unique map of graph
174  A->getColMap(), // column map of A
175  fullblocksize,
176  offset,
177  blockid,
178  nStridedOffset,
179  stridedblocksize) );
180  }
181 
182  // store (un)amalgamation information on current level
183  Set(currentLevel, "UnAmalgamationInfo", amalgamationData);
184  }
185 
186  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
188  AmalgamateMap(const Map& sourceMap, const Matrix& A, RCP<const Map>& amalgamatedMap, Array<LO>& translation) {
189  typedef typename ArrayView<const GO>::size_type size_type;
190  typedef std::map<GO,size_type> container;
191 
192  GO indexBase = sourceMap.getIndexBase();
193  ArrayView<const GO> elementAList = sourceMap.getLocalElementList();
194  size_type numElements = elementAList.size();
195  container filter; // TODO: replace std::set with an object having faster lookup/insert, hashtable for instance
196 
197  GO offset = 0;
198  LO blkSize = A.GetFixedBlockSize() / A.GetStorageBlockSize();
199  if (A.IsView("stridedMaps") == true) {
200  Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
201  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
202  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
203  offset = strMap->getOffset();
204  blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
205  }
206 
207  Array<GO> elementList(numElements);
208  translation.resize(numElements);
209 
210  size_type numRows = 0;
211  for (size_type id = 0; id < numElements; id++) {
212  GO dofID = elementAList[id];
213  GO nodeID = AmalgamationFactory_kokkos::DOFGid2NodeId(dofID, blkSize, offset, indexBase);
214 
215  typename container::iterator it = filter.find(nodeID);
216  if (it == filter.end()) {
217  filter[nodeID] = numRows;
218 
219  translation[id] = numRows;
220  elementList[numRows] = nodeID;
221 
222  numRows++;
223 
224  } else {
225  translation[id] = it->second;
226  }
227  }
228  elementList.resize(numRows);
229 
230  amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
231 
232  }
233 
234  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
236  DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase) {
237  // here, the assumption is, that the node map has the same indexBase as the dof map
238  GlobalOrdinal globalblockid = ((GlobalOrdinal) gid - offset - indexBase) / blockSize + indexBase;
239  return globalblockid;
240  }
241 
242 } //namespace MueLu
243 
244 #endif /* MUELU_AMALGAMATIONFACTORY_KOKKOS_DEF_HPP */
245 
Important warning messages (one line)
Exception indicating invalid cast attempted.
MueLu::DefaultLocalOrdinal LocalOrdinal
void Build(Level &currentLevel) const
Build an object with this factory.
GlobalOrdinal GO
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void DeclareInput(Level &currentLevel) const
Input.
size_type size() const
LocalOrdinal LO
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
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:99
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
minimal container class for storing amalgamation information
void resize(size_type new_size, const value_type &x=value_type())
std::string viewLabel_t
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.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)