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 
55 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
57 #include "MueLu_AmalgamationInfo_kokkos.hpp"
58 
59 namespace MueLu {
60 
61  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
62  RCP<const ParameterList> AmalgamationFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
63  GetValidParameterList() const {
64  RCP<ParameterList> validParamList = rcp(new ParameterList());
65  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
66  return validParamList;
67  }
68 
69  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  void AmalgamationFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
71  DeclareInput(Level &currentLevel) const {
72  Input(currentLevel, "A"); // sub-block from blocked A
73  }
74 
75  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
76  void AmalgamationFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
77  Build(Level &currentLevel) const
78  {
79  FactoryMonitor m(*this, "Build", currentLevel);
80 
81  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
82 
83  LO fullblocksize = 1; // block dim for fixed size blocks
84  GO offset = 0; // global offset of dof gids
85  LO blockid = -1; // block id in strided map
86  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
87  LO stridedblocksize = fullblocksize; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
88  // GO indexBase = A->getRowMap()->getIndexBase(); // index base for maps (unused)
89 
90  // 1) check for blocking/striding information
91 
92  if (A->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
93  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // NOTE: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
94  RCP<const StridedMap> stridedRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
95  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap == Teuchos::null,Exceptions::BadCast,"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
96  fullblocksize = stridedRowMap->getFixedBlockSize();
97  offset = stridedRowMap->getOffset();
98  blockid = stridedRowMap->getStridedBlockId();
99 
100  if (blockid > -1) {
101  std::vector<size_t> stridingInfo = stridedRowMap->getStridingData();
102  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
103  nStridedOffset += stridingInfo[j];
104  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
105 
106  } else {
107  stridedblocksize = fullblocksize;
108  }
109  oldView = A->SwitchToView(oldView);
110  GetOStream(Runtime1) << "AmalagamationFactory::Build():" << " found fullblocksize=" << fullblocksize << " and stridedblocksize=" << stridedblocksize << " from strided maps. offset=" << offset << std::endl;
111 
112  } else {
113  GetOStream(Warnings0) << "AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
114  }
115 
116  // build node row map (uniqueMap) and node column map (nonUniqueMap)
117  // the arrays rowTranslation and colTranslation contain the local node id
118  // given a local dof id. They are only necessary for the CoalesceDropFactory if
119  // fullblocksize > 1
120  RCP<const Map> uniqueMap, nonUniqueMap;
121  RCP<AmalgamationInfo_kokkos> amalgamationData;
122  RCP<Array<LO> > rowTranslation = Teuchos::null;
123  RCP<Array<LO> > colTranslation = Teuchos::null;
124 
125  if (fullblocksize > 1) {
126  // mfh 14 Apr 2015: These need to have different names than
127  // rowTranslation and colTranslation, in order to avoid
128  // shadowing warnings (-Wshadow with GCC). Alternately, it
129  // looks like you could just assign to the existing variables in
130  // this scope, rather than creating new ones.
131  RCP<Array<LO> > theRowTranslation = rcp(new Array<LO>);
132  RCP<Array<LO> > theColTranslation = rcp(new Array<LO>);
133  AmalgamateMap(*(A->getRowMap()), *A, uniqueMap, *theRowTranslation);
134  AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, *theColTranslation);
135 
136  amalgamationData = rcp(new AmalgamationInfo_kokkos(theRowTranslation,
137  theColTranslation,
138  uniqueMap,
139  nonUniqueMap,
140  A->getColMap(),
141  fullblocksize,
142  offset,
143  blockid,
144  nStridedOffset,
145  stridedblocksize) );
146  } else {
147  amalgamationData = rcp(new AmalgamationInfo_kokkos(rowTranslation, // Teuchos::null
148  colTranslation, // Teuchos::null
149  A->getRowMap(), // unique map of graph
150  A->getColMap(), // non-unique map of graph
151  A->getColMap(), // column map of A
152  fullblocksize,
153  offset,
154  blockid,
155  nStridedOffset,
156  stridedblocksize) );
157  }
158 
159  // store (un)amalgamation information on current level
160  Set(currentLevel, "UnAmalgamationInfo", amalgamationData);
161  }
162 
163  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
164  void AmalgamationFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
165  AmalgamateMap(const Map& sourceMap, const Matrix& A, RCP<const Map>& amalgamatedMap, Array<LO>& translation) {
166  typedef typename ArrayView<const GO>::size_type size_type;
167  typedef std::map<GO,size_type> container;
168 
169  GO indexBase = sourceMap.getIndexBase();
170  ArrayView<const GO> elementAList = sourceMap.getNodeElementList();
171  size_type numElements = elementAList.size();
172  container filter; // TODO: replace std::set with an object having faster lookup/insert, hashtable for instance
173 
174  GO offset = 0;
175  LO blkSize = A.GetFixedBlockSize();
176  if (A.IsView("stridedMaps") == true) {
177  Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
178  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
179  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
180  offset = strMap->getOffset();
181  blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
182  }
183 
184  Array<GO> elementList(numElements);
185  translation.resize(numElements);
186 
187  size_type numRows = 0;
188  for (size_type id = 0; id < numElements; id++) {
189  GO dofID = elementAList[id];
190  GO nodeID = AmalgamationFactory_kokkos::DOFGid2NodeId(dofID, blkSize, offset, indexBase);
191 
192  typename container::iterator it = filter.find(nodeID);
193  if (it == filter.end()) {
194  filter[nodeID] = numRows;
195 
196  translation[id] = numRows;
197  elementList[numRows] = nodeID;
198 
199  numRows++;
200 
201  } else {
202  translation[id] = it->second;
203  }
204  }
205  elementList.resize(numRows);
206 
207  amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
208 
209  }
210 
211  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
212  const GlobalOrdinal AmalgamationFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
213  DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase) {
214  // here, the assumption is, that the node map has the same indexBase as the dof map
215  GlobalOrdinal globalblockid = ((GlobalOrdinal) gid - offset - indexBase) / blockSize + indexBase;
216  return globalblockid;
217  }
218 
219 } //namespace MueLu
220 
221 #endif // HAVE_MUELU_KOKKOS_REFACTOR
222 #endif /* MUELU_AMALGAMATIONFACTORY_KOKKOS_DEF_HPP */
223 
Important warning messages (one line)
MueLu::DefaultLocalOrdinal LocalOrdinal
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Description of what is happening (more verbose)