46#ifndef MUELU_AMALGAMATIONFACTORY_DEF_HPP
47#define MUELU_AMALGAMATIONFACTORY_DEF_HPP
49#include <Xpetra_Matrix.hpp>
51#include "MueLu_AmalgamationFactory.hpp"
54#include "MueLu_AmalgamationInfo.hpp"
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;
66 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 Input(currentLevel,
"A");
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
99 LO nStridedOffset = 0;
100 LO stridedblocksize = fullblocksize;
101 LO storageblocksize = A->GetStorageBlockSize();
106 if (A->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
107 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
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();
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]);
121 stridedblocksize = fullblocksize;
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;
129 oldView = A->SwitchToView(oldView);
130 GetOStream(
Runtime1) <<
"AmalagamationFactory::Build():" <<
" found fullblocksize=" << fullblocksize <<
" and stridedblocksize=" << stridedblocksize <<
" from strided maps. offset=" << offset << std::endl;
133 GetOStream(
Warnings0) <<
"AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
141 RCP<const Map> uniqueMap, nonUniqueMap;
142 RCP<AmalgamationInfo> amalgamationData;
143 RCP<Array<LO> > rowTranslation = Teuchos::null;
144 RCP<Array<LO> > colTranslation = Teuchos::null;
146 if (fullblocksize > 1) {
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);
181 Set(currentLevel,
"UnAmalgamationInfo", amalgamationData);
184 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
186 typedef typename ArrayView<const GO>::size_type size_type;
187 typedef std::unordered_map<GO,size_type> container;
189 GO indexBase = sourceMap.getIndexBase();
190 ArrayView<const GO> elementAList = sourceMap.getLocalElementList();
191 size_type numElements = elementAList.size();
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);
200 offset = strMap->getOffset();
201 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
204 Array<GO> elementList(numElements);
205 translation.resize(numElements);
207 size_type numRows = 0;
208 for (size_type
id = 0;
id < numElements;
id++) {
209 GO dofID = elementAList[id];
212 typename container::iterator it = filter.find(nodeID);
213 if (it == filter.end()) {
214 filter[nodeID] = numRows;
216 translation[id] = numRows;
217 elementList[numRows] = nodeID;
222 translation[id] = it->second;
225 elementList.resize(numRows);
227 amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
231 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
235 return globalblockid;
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 ¤tLevel) 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 ¤tLevel) 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.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Runtime1
Description of what is happening (more verbose)