MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_MapTransferFactory_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#ifndef MUELU_MAPTRANSFERFACTORY_DEF_HPP_
48#define MUELU_MAPTRANSFERFACTORY_DEF_HPP_
49
51
52#include <Xpetra_Map.hpp>
53#include <Xpetra_MapFactory.hpp>
54#include <Xpetra_Matrix.hpp>
55
56#include "MueLu_Level.hpp"
58#include "MueLu_Monitor.hpp"
59
60namespace MueLu {
61
62 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
64 RCP<ParameterList> validParamList = rcp(new ParameterList());
65
66 validParamList->setEntry("map: name", Teuchos::ParameterEntry(std::string("")));
67 validParamList->setEntry("map: factory", Teuchos::ParameterEntry(std::string("null")));
68
69 validParamList->set<RCP<const FactoryBase>>("P", Teuchos::null, "Tentative prolongator factory");
70 validParamList->set<std::string>("nullspace vectors: limit to", "all", "Limit the number of nullspace vectors to be used for the map transfer (especially to exclude rotational vectors).");
71
72 return validParamList;
73 }
74
75 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77 const ParameterList & pL = GetParameterList();
78 const std::string mapFactName = pL.get<std::string>("map: factory");
79 const std::string mapName = pL.get<std::string>("map: name");
80
81 if (fineLevel.GetLevelID() == 0)
82 {
83 // Not needed, if the map is provided as user data
84 fineLevel.DeclareInput(mapName, NoFactory::get(), this);
85 }
86 else
87 {
88 // check whether user has provided a specific name for the MapFactory
89 if (mapFactName == "" || mapFactName == "NoFactory")
90 mapFact_ = MueLu::NoFactory::getRCP();
91 else if (mapFactName != "null")
92 mapFact_ = coarseLevel.GetFactoryManager()->GetFactory(mapFactName);
93
94 // request map generated by mapFact_
95 fineLevel.DeclareInput(mapName, mapFact_.get(), this);
96 }
97
98 // request Ptent
99 // note that "P" provided by the user (through XML file) is supposed to be of type TentativePFactory
100 Teuchos::RCP<const FactoryBase> tentPFact = GetFactory("P");
101 if (tentPFact == Teuchos::null)
102 tentPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
103 coarseLevel.DeclareInput("P", tentPFact.get(), this);
104 }
105
106 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
108 Monitor m(*this, "MapTransferFactory");
109
110 const ParameterList & pL = GetParameterList();
111 const std::string mapName = pL.get<std::string>("map: name");
112 const int maxNumProlongCols = GetLimitOfProlongatorColumns(pL);
113
114 // fetch map from level
115 RCP<const Map> transferMap = Teuchos::null;
116 if (fineLevel.GetLevelID() == 0) {
117 transferMap = fineLevel.Get<RCP<const Map>>(mapName, NoFactory::get());
118 } else {
119 if (fineLevel.IsAvailable(mapName, mapFact_.get()) == false)
120 GetOStream(Runtime0) << "MapTransferFactory::Build: User provided map \"" << mapName << "\" not found in Level class on level " << fineLevel.GetLevelID() << "." << std::endl;
121 transferMap = fineLevel.Get<RCP<const Map>>(mapName, mapFact_.get());
122 }
123
124 // Get default tentative prolongator factory
125 // Getting it that way ensures that the same factory instance will be used for both SaPFactory and NullspaceFactory.
126 RCP<const FactoryBase> tentPFact = GetFactory("P");
127 if (tentPFact == Teuchos::null)
128 tentPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
129 TEUCHOS_TEST_FOR_EXCEPTION(!coarseLevel.IsAvailable("P", tentPFact.get()), Exceptions::RuntimeError,
130 "MueLu::MapTransferFactory::Build(): P (generated by TentativePFactory) not available.");
131 RCP<Matrix> Ptent = coarseLevel.Get<RCP<Matrix> >("P", tentPFact.get());
132
133 // loop over local rows of Ptent and figure out the corresponding coarse GIDs
134 Array<GO > coarseMapGids;
135 RCP<const Map> prolongColMap = Ptent->getColMap();
136 GO gRowID = -1;
137 int numColEntries = 0;
138 for (size_t row = 0; row < Ptent->getLocalNumRows(); ++row) {
139 gRowID = Ptent->getRowMap()->getGlobalElement(row);
140
141 if (transferMap->isNodeGlobalElement(gRowID)) {
142 Teuchos::ArrayView<const LO> indices;
143 Teuchos::ArrayView<const SC> vals;
144 Ptent->getLocalRowView(row, indices, vals);
145
146 numColEntries = as<int>(indices.size());
147 if (maxNumProlongCols > 0)
148 numColEntries = std::min(numColEntries, maxNumProlongCols);
149
150 for (size_t col = 0; col < as<size_t>(numColEntries); ++col) {
151 // mark all (selected) columns in Ptent(gRowID,*) to be coarse Dofs of next level transferMap
152 GO gcid = prolongColMap->getGlobalElement(indices[col]);
153 coarseMapGids.push_back(gcid);
154 }
155 }
156 }
157
158 // build coarse version of the input map
159 const GO INVALID = Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
160 std::sort(coarseMapGids.begin(), coarseMapGids.end());
161 coarseMapGids.erase(std::unique(coarseMapGids.begin(), coarseMapGids.end()), coarseMapGids.end());
162 RCP<const Map> coarseTransferMap = MapFactory::Build(prolongColMap->lib(), INVALID, coarseMapGids(),
163 prolongColMap->getIndexBase(), prolongColMap->getComm());
164
165 // store map in coarse level
166 if (fineLevel.GetLevelID() == 0)
167 {
168 const std::string mapFactName = pL.get<std::string>("map: factory");
169 RCP<const FactoryBase> mapFact = coarseLevel.GetFactoryManager()->GetFactory(mapFactName);
170 coarseLevel.Set(mapName, coarseTransferMap, mapFact.get());
171 }
172 else
173 coarseLevel.Set(mapName, coarseTransferMap, mapFact_.get());
174
175 }
176
177 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
179 {
180 const std::string useTheseNspVectors = pL.get<std::string>("nullspace vectors: limit to");
181
182 // Leave right away, if no limit is prescribed by the user
183 if (useTheseNspVectors == "all" || useTheseNspVectors == "")
184 return -1;
185
186 // Simplify? Maybe replace by boolean flag "nullspace: exclude rotations"
187 int maxNumProlongCols = -1;
188 if (useTheseNspVectors == "translations")
189 maxNumProlongCols = 1;
190 else
191 TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::InvalidArgument, "Unknown subset of nullspace vectors to be used, when performing a map transfer.")
192
193 return maxNumProlongCols;
194 }
195
196} // namespace MueLu
197
198#endif /* MUELU_MAPTRANSFERFACTORY_DEF_HPP_ */
Exception throws to report invalid user entry.
Exception throws to report errors in the internal logical of the program.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
int GetLimitOfProlongatorColumns(const ParameterList &pL) const
Get the max number of entries per row of P to be considered for map transfer.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const override
Input.
void Build(Level &fineLevel, Level &coarseLevel) const override
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const override
Input.
Timer to be used in non-factories.
static const RCP< const NoFactory > getRCP()
Static Get() functions.
static const NoFactory * get()
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.