MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_AmalgamationInfo_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_DEF_HPP_
54#define MUELU_AMALGAMATIONINFO_DEF_HPP_
55
56#include <Xpetra_MapFactory.hpp>
57#include <Xpetra_Vector.hpp>
58
60#include "MueLu_Exceptions.hpp"
61#include "MueLu_Aggregates.hpp"
62
63namespace MueLu {
64
65 template <class LocalOrdinal, class GlobalOrdinal, class Node>
67 Teuchos::ArrayRCP<LocalOrdinal>& aggStart, Teuchos::ArrayRCP<GlobalOrdinal>& aggToRowMap) const {
68 const int myPid = aggregates.GetMap()->getComm()->getRank();
69 Teuchos::ArrayView<const GO> nodeGlobalElts = aggregates.GetMap()->getLocalElementList();
70 Teuchos::ArrayRCP<LO> procWinner = aggregates.GetProcWinner()->getDataNonConst(0);
71 Teuchos::ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
72 const LO size = procWinner.size();
73 const GO numAggregates = aggregates.GetNumAggregates();
74
75 Array<LO> sizes(numAggregates);
76 if (stridedblocksize_ == 1) {
77 for (LO lnode = 0; lnode < size; ++lnode) {
78 LO myAgg = vertex2AggId[lnode];
79 if (procWinner[lnode] == myPid)
80 sizes[myAgg] += 1;
81 }
82 } else {
83 for (LO lnode = 0; lnode < size; ++lnode) {
84 LO myAgg = vertex2AggId[lnode];
85 if (procWinner[lnode] == myPid) {
86 GO gnodeid = nodeGlobalElts[lnode];
87 for (LocalOrdinal k = 0; k < stridedblocksize_; k++) {
88 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gnodeid,k);
89 if (columnMap_->isNodeGlobalElement(gDofIndex))
90 sizes[myAgg] += 1;
91 }
92 }
93 }
94 }
95 aggStart = ArrayRCP<LO>(numAggregates+1,0);
96 aggStart[0] = Teuchos::ScalarTraits<LO>::zero();
97 for (GO i=0; i<numAggregates; ++i) {
98 aggStart[i+1] = aggStart[i] + sizes[i];
99 }
100 aggToRowMap = ArrayRCP<GO>(aggStart[numAggregates],0);
101
102 // count, how many dofs have been recorded for each aggregate so far
103 Array<LO> numDofs(numAggregates, 0); // empty array with number of Dofs for each aggregate
104
105 if (stridedblocksize_ == 1) {
106 for (LO lnode = 0; lnode < size; ++lnode) {
107 LO myAgg = vertex2AggId[lnode];
108 if (procWinner[lnode] == myPid) {
109 aggToRowMap[ aggStart[myAgg] + numDofs[myAgg] ] = ComputeGlobalDOF(nodeGlobalElts[lnode]);
110 ++(numDofs[myAgg]);
111 }
112 }
113 } else {
114 for (LO lnode = 0; lnode < size; ++lnode) {
115 LO myAgg = vertex2AggId[lnode];
116
117 if (procWinner[lnode] == myPid) {
118 GO gnodeid = nodeGlobalElts[lnode];
119 for (LocalOrdinal k = 0; k < stridedblocksize_; k++) {
120 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gnodeid,k);
121 if (columnMap_->isNodeGlobalElement(gDofIndex)) {
122 aggToRowMap[ aggStart[myAgg] + numDofs[myAgg] ] = gDofIndex;
123 ++(numDofs[myAgg]);
124 }
125 }
126 }
127 }
128 }
129 // todo plausibility check: entry numDofs[k] == aggToRowMap[k].size()
130
131 } //UnamalgamateAggregates
132
133 template <class LocalOrdinal, class GlobalOrdinal, class Node>
135 Teuchos::ArrayRCP<LO>& aggStart, Teuchos::ArrayRCP<LO>& aggToRowMap) const {
136
137 int myPid = aggregates.GetMap()->getComm()->getRank();
138 Teuchos::ArrayView<const GO> nodeGlobalElts = aggregates.GetMap()->getLocalElementList();
139
140 Teuchos::ArrayRCP<LO> procWinner = aggregates.GetProcWinner() ->getDataNonConst(0);
141 Teuchos::ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
142 const GO numAggregates = aggregates.GetNumAggregates();
143
144
145 // FIXME: Do we need to compute size here? Or can we use existing?
146 const LO size = procWinner.size();
147
148 std::vector<LO> sizes(numAggregates);
149 if (stridedblocksize_ == 1) {
150 for (LO lnode = 0; lnode < size; lnode++)
151 if (procWinner[lnode] == myPid)
152 sizes[vertex2AggId[lnode]]++;
153 } else {
154 for (LO lnode = 0; lnode < size; lnode++)
155 if (procWinner[lnode] == myPid) {
156 GO nodeGID = nodeGlobalElts[lnode];
157
158 for (LO k = 0; k < stridedblocksize_; k++) {
159 GO GID = ComputeGlobalDOF(nodeGID, k);
160 if (columnMap_->isNodeGlobalElement(GID))
161 sizes[vertex2AggId[lnode]]++;
162 }
163 }
164 }
165
166 aggStart = ArrayRCP<LO>(numAggregates+1); // FIXME: useless initialization with zeros
167 aggStart[0] = 0;
168 for (GO i = 0; i < numAggregates; i++)
169 aggStart[i+1] = aggStart[i] + sizes[i];
170
171 aggToRowMap = ArrayRCP<LO>(aggStart[numAggregates], 0);
172
173 // count, how many dofs have been recorded for each aggregate so far
174 Array<LO> numDofs(numAggregates, 0); // empty array with number of DOFs for each aggregate
175 if (stridedblocksize_ == 1) {
176 for (LO lnode = 0; lnode < size; ++lnode)
177 if (procWinner[lnode] == myPid) {
178 LO myAgg = vertex2AggId[lnode];
179 aggToRowMap[aggStart[myAgg] + numDofs[myAgg]] = lnode;
180 numDofs[myAgg]++;
181 }
182 } else {
183 for (LO lnode = 0; lnode < size; ++lnode)
184 if (procWinner[lnode] == myPid) {
185 LO myAgg = vertex2AggId[lnode];
186 GO nodeGID = nodeGlobalElts[lnode];
187
188 for (LO k = 0; k < stridedblocksize_; k++) {
189 GO GID = ComputeGlobalDOF(nodeGID, k);
190 if (columnMap_->isNodeGlobalElement(GID)) {
191 aggToRowMap[aggStart[myAgg] + numDofs[myAgg]] = lnode*stridedblocksize_ + k;
192 numDofs[myAgg]++;
193 }
194 }
195 }
196 }
197 // todo plausibility check: entry numDofs[k] == aggToRowMap[k].size()
198
199 } //UnamalgamateAggregatesLO
200
201 template <class LocalOrdinal, class GlobalOrdinal, class Node>
203 const VerbLevel verbLevel) const
204 {
205 if (!(verbLevel & Debug))
206 return;
207
208 out << "AmalgamationInfo -- Striding information:"
209 << "\n fullBlockSize = " << fullblocksize_
210 << "\n blockID = " << blockid_
211 << "\n stridingOffset = " << nStridedOffset_
212 << "\n stridedBlockSize = " << stridedblocksize_
213 << "\n indexBase = " << indexBase_
214 << std::endl;
215
216 out << "AmalgamationInfo -- DOFs to nodes mapping:\n"
217 << " Mapping of row DOFs to row nodes:" << *rowTranslation_()
218 << "\n\n Mapping of column DOFs to column nodes:" << *colTranslation_()
219 << std::endl;
220
221 out << "AmalgamationInfo -- row node map:" << std::endl;
222 nodeRowMap_->describe(out, Teuchos::VERB_EXTREME);
223
224 out << "AmalgamationInfo -- column node map:" << std::endl;
225 nodeColMap_->describe(out, Teuchos::VERB_EXTREME);
226 }
227
229
230 template <class LocalOrdinal, class GlobalOrdinal, class Node>
231 RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>::ComputeUnamalgamatedImportDofMap(const Aggregates& aggregates) const {
232 Teuchos::RCP<const Map> nodeMap = aggregates.GetMap();
233
234 Teuchos::RCP<std::vector<GO> > myDofGids = Teuchos::rcp(new std::vector<GO>);
235 Teuchos::ArrayView<const GO> gEltList = nodeMap->getLocalElementList();
236 LO nodeElements = Teuchos::as<LO>(nodeMap->getLocalNumElements());
237 if (stridedblocksize_ == 1) {
238 for (LO n = 0; n<nodeElements; n++) {
239 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gEltList[n]);
240 myDofGids->push_back(gDofIndex);
241 }
242 } else {
243 for (LO n = 0; n<nodeElements; n++) {
244 for (LocalOrdinal k = 0; k < stridedblocksize_; k++) {
245 GlobalOrdinal gDofIndex = ComputeGlobalDOF(gEltList[n],k);
246 if (columnMap_->isNodeGlobalElement(gDofIndex))
247 myDofGids->push_back(gDofIndex);
248 }
249 }
250 }
251
252 Teuchos::ArrayRCP<GO> arr_myDofGids = Teuchos::arcp( myDofGids );
253 Teuchos::RCP<Map> importDofMap = MapFactory::Build(aggregates.GetMap()->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), arr_myDofGids(), aggregates.GetMap()->getIndexBase(), aggregates.GetMap()->getComm());
254 return importDofMap;
255 }
256
258
259 template <class LocalOrdinal, class GlobalOrdinal, class Node>
261 GlobalOrdinal gDofIndex = offset_ + (gNodeID-indexBase_)*fullblocksize_ + nStridedOffset_ + k + indexBase_;
262 return gDofIndex;
263 }
264
265 template <class LocalOrdinal, class GlobalOrdinal, class Node>
267 LocalOrdinal lDofIndex = lNodeID*fullblocksize_ + k;
268 return lDofIndex;
269 }
270
271
272 template <class LocalOrdinal, class GlobalOrdinal, class Node>
274 return (ldofID - ldofID%fullblocksize_) / fullblocksize_;
275 }
276
277} //namespace
278
279
280#endif /* MUELU_AMALGAMATIONINFO_DEF_HPP_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Container class for aggregation information.
const RCP< const Map > GetMap() const
returns (overlapping) map of aggregate/node distribution
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
LO GetNumAggregates() const
returns the number of aggregates of the current processor. Note: could/should be renamed to GetNumLoc...
LO ComputeLocalDOF(LocalOrdinal const &lNodeID, LocalOrdinal const &k) const
ComputeLocalDOF return locbal dof id associated with local node id lNodeID and dof index k.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void UnamalgamateAggregatesLO(const Aggregates &aggregates, Teuchos::ArrayRCP< LocalOrdinal > &aggStart, Teuchos::ArrayRCP< LO > &aggToRowMap) const
GO ComputeGlobalDOF(GO const &gNodeID, LO const &k=0) const
ComputeGlobalDOF.
Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ComputeUnamalgamatedImportDofMap(const Aggregates &aggregates) const
ComputeUnamalgamatedImportDofMap build overlapping dof row map from aggregates needed for overlapping...
void UnamalgamateAggregates(const Aggregates &aggregates, Teuchos::ArrayRCP< LocalOrdinal > &aggStart, Teuchos::ArrayRCP< GlobalOrdinal > &aggToRowMap) const
UnamalgamateAggregates.
LO ComputeLocalNode(LocalOrdinal const &ldofID) const
Namespace for MueLu classes and methods.
@ Debug
Print additional debugging information.