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