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");
81 LO nStridedOffset = 0;
82 LO stridedblocksize = fullblocksize;
87 if (A->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
88 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
89 RCP<const StridedMap> stridedRowMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap());
90 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
91 fullblocksize = stridedRowMap->getFixedBlockSize();
92 offset = stridedRowMap->getOffset();
93 blockid = stridedRowMap->getStridedBlockId();
96 std::vector<size_t> stridingInfo = stridedRowMap->getStridingData();
97 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
98 nStridedOffset += stridingInfo[j];
99 stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
102 stridedblocksize = fullblocksize;
104 oldView = A->SwitchToView(oldView);
105 GetOStream(
Runtime1) <<
"AmalagamationFactory::Build():" <<
" found fullblocksize=" << fullblocksize <<
" and stridedblocksize=" << stridedblocksize <<
" from strided maps. offset=" << offset << std::endl;
108 GetOStream(
Warnings0) <<
"AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
115 RCP<const Map> uniqueMap, nonUniqueMap;
116 RCP<AmalgamationInfo> amalgamationData;
117 RCP<Array<LO> > rowTranslation = Teuchos::null;
118 RCP<Array<LO> > colTranslation = Teuchos::null;
120 if (fullblocksize > 1) {
126 RCP<Array<LO> > theRowTranslation = rcp(
new Array<LO>);
127 RCP<Array<LO> > theColTranslation = rcp(
new Array<LO>);
128 AmalgamateMap(*(A->getRowMap()), *A, uniqueMap, *theRowTranslation);
129 AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, *theColTranslation);
155 Set(currentLevel,
"UnAmalgamationInfo", amalgamationData);
158 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
160 typedef typename ArrayView<const GO>::size_type size_type;
161 typedef std::map<GO,size_type> container;
163 GO indexBase = sourceMap.getIndexBase();
164 ArrayView<const GO> elementAList = sourceMap.getNodeElementList();
165 size_type numElements = elementAList.size();
169 LO blkSize = A.GetFixedBlockSize();
170 if (A.IsView(
"stridedMaps") ==
true) {
171 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
172 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
174 offset = strMap->getOffset();
175 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
178 Array<GO> elementList(numElements);
179 translation.resize(numElements);
181 size_type numRows = 0;
182 for (size_type
id = 0;
id < numElements;
id++) {
183 GO dofID = elementAList[id];
186 typename container::iterator it = filter.find(nodeID);
187 if (it == filter.end()) {
188 filter[nodeID] = numRows;
190 translation[id] = numRows;
191 elementList[numRows] = nodeID;
196 translation[id] = it->second;
199 elementList.resize(numRows);
201 amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
205 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
208 GlobalOrdinal globalblockid = ((GlobalOrdinal) gid - offset - indexBase) / blockSize + indexBase;
209 return globalblockid;
Important warning messages (one line)
Exception indicating invalid cast attempted.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
Class that holds all level-specific information.
void DeclareInput(Level ¤tLevel) const
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.
void Build(Level ¤tLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
minimal container class for storing amalgamation information