MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_AmalgamationInfo_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/*
47 * MueLu_AmalgamationInfo_def.hpp
48 *
49 * Created on: Mar 28, 2012
50 * Author: wiesner
51 */
52
53#ifndef MUELU_AMALGAMATIONINFO_KOKKOS_DEF_HPP_
54#define MUELU_AMALGAMATIONINFO_KOKKOS_DEF_HPP_
55
56#include <Xpetra_MapFactory.hpp>
57#include <Xpetra_Vector.hpp>
58
59#include "MueLu_Exceptions.hpp"
60
62#include "MueLu_Aggregates_kokkos.hpp"
63
64namespace MueLu {
65
66 template <class LocalOrdinal, class GlobalOrdinal, class Node>
69 Teuchos::ArrayRCP<LocalOrdinal>& aggStart,
70 Teuchos::ArrayRCP<GlobalOrdinal>& aggToRowMap) const {
71
72 int myPid = aggregates.GetMap()->getComm()->getRank();
73 Teuchos::ArrayView<const GO> nodeGlobalElts = aggregates.GetMap()->getLocalElementList();
74 Teuchos::ArrayRCP<LO> procWinner = aggregates.GetProcWinner()->getDataNonConst(0);
75 Teuchos::ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
76 LO size = procWinner.size();
77 GO numAggregates = aggregates.GetNumAggregates();
78
79 std::vector<LO> sizes(numAggregates);
80 if (stridedblocksize_ == 1) {
81 for (LO lnode = 0; lnode < size; ++lnode) {
82 LO myAgg = vertex2AggId[lnode];
83 if (procWinner[lnode] == myPid)
84 sizes[myAgg] += 1;
85 }
86 } else {
87 for (LO lnode = 0; lnode < size; ++lnode) {
88 LO myAgg = vertex2AggId[lnode];
89 if (procWinner[lnode] == myPid) {
90 GO gnodeid = nodeGlobalElts[lnode];
91 for (LocalOrdinal k = 0; k < stridedblocksize_; k++) {
92 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gnodeid,k);
93 if (columnMap_->isNodeGlobalElement(gDofIndex))
94 sizes[myAgg] += 1;
95 }
96 }
97 }
98 }
99 aggStart = ArrayRCP<LO>(numAggregates+1,0);
100 aggStart[0]=0;
101 for (GO i=0; i<numAggregates; ++i) {
102 aggStart[i+1] = aggStart[i] + sizes[i];
103 }
104 aggToRowMap = ArrayRCP<GO>(aggStart[numAggregates],0);
105
106 // count, how many dofs have been recorded for each aggregate so far
107 Array<LO> numDofs(numAggregates, 0); // empty array with number of Dofs for each aggregate
108
109 if (stridedblocksize_ == 1) {
110 for (LO lnode = 0; lnode < size; ++lnode) {
111 LO myAgg = vertex2AggId[lnode];
112 if (procWinner[lnode] == myPid) {
113 aggToRowMap[ aggStart[myAgg] + numDofs[myAgg] ] = ComputeGlobalDOF(nodeGlobalElts[lnode]);
114 ++(numDofs[myAgg]);
115 }
116 }
117 } else {
118 for (LO lnode = 0; lnode < size; ++lnode) {
119 LO myAgg = vertex2AggId[lnode];
120
121 if (procWinner[lnode] == myPid) {
122 GO gnodeid = nodeGlobalElts[lnode];
123 for (LocalOrdinal k = 0; k < stridedblocksize_; k++) {
124 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gnodeid,k);
125 if (columnMap_->isNodeGlobalElement(gDofIndex)) {
126 aggToRowMap[ aggStart[myAgg] + numDofs[myAgg] ] = gDofIndex;
127 ++(numDofs[myAgg]);
128 }
129 }
130 }
131 }
132 }
133 // todo plausibility check: entry numDofs[k] == aggToRowMap[k].size()
134
135 } //UnamalgamateAggregates
136
137 template <class LocalOrdinal, class GlobalOrdinal, class Node>
140 Teuchos::ArrayRCP<LO>& aggStart,
141 Teuchos::ArrayRCP<LO>& aggToRowMap) const {
142
143 int myPid = aggregates.GetMap()->getComm()->getRank();
144 Teuchos::ArrayView<const GO> nodeGlobalElts = aggregates.GetMap()->getLocalElementList();
145
146 Teuchos::ArrayRCP<LO> procWinner = aggregates.GetProcWinner() ->getDataNonConst(0);
147 Teuchos::ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
148 const GO numAggregates = aggregates.GetNumAggregates();
149
150
151 // FIXME: Do we need to compute size here? Or can we use existing?
152 LO size = procWinner.size();
153
154 std::vector<LO> sizes(numAggregates);
155 if (stridedblocksize_ == 1) {
156 for (LO lnode = 0; lnode < size; lnode++)
157 if (procWinner[lnode] == myPid)
158 sizes[vertex2AggId[lnode]]++;
159 } else {
160 for (LO lnode = 0; lnode < size; lnode++)
161 if (procWinner[lnode] == myPid) {
162 GO nodeGID = nodeGlobalElts[lnode];
163
164 for (LO k = 0; k < stridedblocksize_; k++) {
165 GO GID = ComputeGlobalDOF(nodeGID, k);
166 if (columnMap_->isNodeGlobalElement(GID))
167 sizes[vertex2AggId[lnode]]++;
168 }
169 }
170 }
171
172 aggStart = ArrayRCP<LO>(numAggregates+1); // FIXME: useless initialization with zeros
173 aggStart[0] = 0;
174 for (GO i = 0; i < numAggregates; i++)
175 aggStart[i+1] = aggStart[i] + sizes[i];
176
177 aggToRowMap = ArrayRCP<LO>(aggStart[numAggregates], 0);
178
179 // count, how many dofs have been recorded for each aggregate so far
180 Array<LO> numDofs(numAggregates, 0); // empty array with number of DOFs for each aggregate
181 if (stridedblocksize_ == 1) {
182 for (LO lnode = 0; lnode < size; ++lnode)
183 if (procWinner[lnode] == myPid) {
184 LO myAgg = vertex2AggId[lnode];
185 aggToRowMap[aggStart[myAgg] + numDofs[myAgg]] = lnode;
186 numDofs[myAgg]++;
187 }
188 } else {
189 for (LO lnode = 0; lnode < size; ++lnode)
190 if (procWinner[lnode] == myPid) {
191 LO myAgg = vertex2AggId[lnode];
192 GO nodeGID = nodeGlobalElts[lnode];
193
194 for (LO k = 0; k < stridedblocksize_; k++) {
195 GO GID = ComputeGlobalDOF(nodeGID, k);
196 if (columnMap_->isNodeGlobalElement(GID)) {
197 aggToRowMap[aggStart[myAgg] + numDofs[myAgg]] = lnode*stridedblocksize_ + k;
198 numDofs[myAgg]++;
199 }
200 }
201 }
202 }
203 // todo plausibility check: entry numDofs[k] == aggToRowMap[k].size()
204
205 } //UnamalgamateAggregates
206
208
209 template <class LocalOrdinal, class GlobalOrdinal, class Node>
210 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > AmalgamationInfo_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
212
213 Teuchos::RCP<const Map> nodeMap = aggregates.GetMap();
214
215 Teuchos::RCP<std::vector<GO> > myDofGids = Teuchos::rcp(new std::vector<GO>);
216 Teuchos::ArrayView<const GO> gEltList = nodeMap->getLocalElementList();
217 LO nodeElements = Teuchos::as<LO>(nodeMap->getLocalNumElements());
218 if (stridedblocksize_ == 1) {
219 for (LO n = 0; n<nodeElements; n++) {
220 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gEltList[n]);
221 myDofGids->push_back(gDofIndex);
222 }
223 } else {
224 for (LO n = 0; n<nodeElements; n++) {
225 for (LocalOrdinal k = 0; k < stridedblocksize_; k++) {
226 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gEltList[n],k);
227 if (columnMap_->isNodeGlobalElement(gDofIndex))
228 myDofGids->push_back(gDofIndex);
229 }
230 }
231 }
232
233 Teuchos::ArrayRCP<GO> arr_myDofGids = Teuchos::arcp( myDofGids );
234 Teuchos::RCP<Map> importDofMap = MapFactory::Build(aggregates.GetMap()->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), arr_myDofGids(), aggregates.GetMap()->getIndexBase(), aggregates.GetMap()->getComm());
235 return importDofMap;
236 }
237
239
240 template <class LocalOrdinal, class GlobalOrdinal, class Node>
242 ComputeGlobalDOF(GlobalOrdinal const &gNodeID, LocalOrdinal const &k) const {
243 // here, the assumption is, that the node map has the same indexBase as the dof map
244 // this is the node map index base this is the dof map index base
245 GlobalOrdinal gDofIndex = offset_ + (gNodeID-indexBase_)*fullblocksize_ + nStridedOffset_ + k + indexBase_;
246 return gDofIndex;
247 }
248
249} //namespace
250
251
252#endif /* MUELU_AMALGAMATIONINFO_KOKKOS_DEF_HPP_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
void UnamalgamateAggregatesLO(const Aggregates_kokkos &aggregates, Teuchos::ArrayRCP< LocalOrdinal > &aggStart, Teuchos::ArrayRCP< LO > &aggToRowMap) const
void UnamalgamateAggregates(const Aggregates_kokkos &aggregates, Teuchos::ArrayRCP< LocalOrdinal > &aggStart, Teuchos::ArrayRCP< GlobalOrdinal > &aggToRowMap) const
UnamalgamateAggregates.
GO ComputeGlobalDOF(GO const &gNodeID, LO const &k=0) const
ComputeGlobalDOF return global dof id associated with global node id gNodeID and dof index k.
Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ComputeUnamalgamatedImportDofMap(const Aggregates_kokkos &aggregates) const
ComputeUnamalgamatedImportDofMap build overlapping dof row map from aggregates needed for overlapping...
Namespace for MueLu classes and methods.