MueLu Version of the Day
Loading...
Searching...
No Matches
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
58namespace 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
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
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 DeclareInput(Level &currentLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level &currentLevel) const
Build an object with this factory.
minimal container class for storing amalgamation information
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Runtime1
Description of what is happening (more verbose)