46#ifndef MUELU_AMALGAMATIONFACTORY_KOKKOS_DEF_HPP
47#define MUELU_AMALGAMATIONFACTORY_KOKKOS_DEF_HPP
49#include <Xpetra_Matrix.hpp>
56#include "MueLu_AmalgamationInfo_kokkos.hpp"
60 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
63 RCP<ParameterList> validParamList = rcp(
new ParameterList());
64 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
65 return validParamList;
68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 Input(currentLevel,
"A");
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
102 LO nStridedOffset = 0;
103 LO stridedblocksize = fullblocksize;
104 LO storageblocksize = A->GetStorageBlockSize();
109 if (A->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
110 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
111 RCP<const StridedMap> stridedRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
112 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
113 fullblocksize = stridedRowMap->getFixedBlockSize();
114 offset = stridedRowMap->getOffset();
115 blockid = stridedRowMap->getStridedBlockId();
118 std::vector<size_t> stridingInfo = stridedRowMap->getStridingData();
119 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
120 nStridedOffset += stridingInfo[j];
121 stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
124 stridedblocksize = fullblocksize;
128 TEUCHOS_TEST_FOR_EXCEPTION(fullblocksize % storageblocksize != 0,
Exceptions::RuntimeError,
"AmalgamationFactory::Build(): fullblocksize needs to be a multiple of A->GetStorageBlockSize()");
129 fullblocksize /= storageblocksize;
130 stridedblocksize /= storageblocksize;
132 oldView = A->SwitchToView(oldView);
133 GetOStream(
Runtime1) <<
"AmalagamationFactory::Build():" <<
" found fullblocksize=" << fullblocksize <<
" and stridedblocksize=" << stridedblocksize <<
" from strided maps. offset=" << offset << std::endl;
136 GetOStream(
Warnings0) <<
"AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
143 RCP<const Map> uniqueMap, nonUniqueMap;
144 RCP<AmalgamationInfo_kokkos> amalgamationData;
145 RCP<Array<LO> > rowTranslation = Teuchos::null;
146 RCP<Array<LO> > colTranslation = Teuchos::null;
148 if (fullblocksize > 1) {
154 RCP<Array<LO> > theRowTranslation = rcp(
new Array<LO>);
155 RCP<Array<LO> > theColTranslation = rcp(
new Array<LO>);
156 AmalgamateMap(*(A->getRowMap()), *A, uniqueMap, *theRowTranslation);
157 AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, *theColTranslation);
183 Set(currentLevel,
"UnAmalgamationInfo", amalgamationData);
186 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
188 AmalgamateMap(
const Map& sourceMap,
const Matrix& A, RCP<const Map>& amalgamatedMap, Array<LO>& translation) {
189 typedef typename ArrayView<const GO>::size_type size_type;
190 typedef std::map<GO,size_type> container;
192 GO indexBase = sourceMap.getIndexBase();
193 ArrayView<const GO> elementAList = sourceMap.getLocalElementList();
194 size_type numElements = elementAList.size();
198 LO blkSize = A.GetFixedBlockSize() / A.GetStorageBlockSize();
199 if (A.IsView(
"stridedMaps") ==
true) {
200 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
201 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
203 offset = strMap->getOffset();
204 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
207 Array<GO> elementList(numElements);
208 translation.resize(numElements);
210 size_type numRows = 0;
211 for (size_type
id = 0;
id < numElements;
id++) {
212 GO dofID = elementAList[id];
215 typename container::iterator it = filter.find(nodeID);
216 if (it == filter.end()) {
217 filter[nodeID] = numRows;
219 translation[id] = numRows;
220 elementList[numRows] = nodeID;
225 translation[id] = it->second;
228 elementList.resize(numRows);
230 amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
234 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
239 return globalblockid;
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
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 DeclareInput(Level ¤tLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level ¤tLevel) const
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)