MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
60 namespace 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  // Not needed, if the map is provided as user data
83  fineLevel.DeclareInput(mapName, NoFactory::get(), this);
84  } else {
85  // check whether user has provided a specific name for the MapFactory
86  if (mapFactName == "" || mapFactName == "NoFactory")
87  mapFact_ = MueLu::NoFactory::getRCP();
88  else if (mapFactName != "null")
89  mapFact_ = coarseLevel.GetFactoryManager()->GetFactory(mapFactName);
90 
91  // request map generated by mapFact_
92  fineLevel.DeclareInput(mapName, mapFact_.get(), this);
93  }
94 
95  // request Ptent
96  // note that "P" provided by the user (through XML file) is supposed to be of type TentativePFactory
97  Teuchos::RCP<const FactoryBase> tentPFact = GetFactory("P");
98  if (tentPFact == Teuchos::null)
99  tentPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
100  coarseLevel.DeclareInput("P", tentPFact.get(), this);
101 }
102 
103 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
105  Monitor m(*this, "MapTransferFactory");
106 
107  const ParameterList& pL = GetParameterList();
108  const std::string mapName = pL.get<std::string>("map: name");
109  const int maxNumProlongCols = GetLimitOfProlongatorColumns(pL);
110 
111  // fetch map from level
112  RCP<const Map> transferMap = Teuchos::null;
113  if (fineLevel.GetLevelID() == 0) {
114  transferMap = fineLevel.Get<RCP<const Map>>(mapName, NoFactory::get());
115  } else {
116  if (fineLevel.IsAvailable(mapName, mapFact_.get()) == false)
117  GetOStream(Runtime0) << "MapTransferFactory::Build: User provided map \"" << mapName << "\" not found in Level class on level " << fineLevel.GetLevelID() << "." << std::endl;
118  transferMap = fineLevel.Get<RCP<const Map>>(mapName, mapFact_.get());
119  }
120 
121  // Get default tentative prolongator factory
122  // Getting it that way ensures that the same factory instance will be used for both SaPFactory and NullspaceFactory.
123  RCP<const FactoryBase> tentPFact = GetFactory("P");
124  if (tentPFact == Teuchos::null)
125  tentPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
126  TEUCHOS_TEST_FOR_EXCEPTION(!coarseLevel.IsAvailable("P", tentPFact.get()), Exceptions::RuntimeError,
127  "MueLu::MapTransferFactory::Build(): P (generated by TentativePFactory) not available.");
128  RCP<Matrix> Ptent = coarseLevel.Get<RCP<Matrix>>("P", tentPFact.get());
129 
130  // loop over local rows of Ptent and figure out the corresponding coarse GIDs
131  Array<GO> coarseMapGids;
132  RCP<const Map> prolongColMap = Ptent->getColMap();
133  GO gRowID = -1;
134  int numColEntries = 0;
135  for (size_t row = 0; row < Ptent->getLocalNumRows(); ++row) {
136  gRowID = Ptent->getRowMap()->getGlobalElement(row);
137 
138  if (transferMap->isNodeGlobalElement(gRowID)) {
141  Ptent->getLocalRowView(row, indices, vals);
142 
143  numColEntries = as<int>(indices.size());
144  if (maxNumProlongCols > 0)
145  numColEntries = std::min(numColEntries, maxNumProlongCols);
146 
147  for (size_t col = 0; col < as<size_t>(numColEntries); ++col) {
148  // mark all (selected) columns in Ptent(gRowID,*) to be coarse Dofs of next level transferMap
149  GO gcid = prolongColMap->getGlobalElement(indices[col]);
150  coarseMapGids.push_back(gcid);
151  }
152  }
153  }
154 
155  // build coarse version of the input map
157  std::sort(coarseMapGids.begin(), coarseMapGids.end());
158  coarseMapGids.erase(std::unique(coarseMapGids.begin(), coarseMapGids.end()), coarseMapGids.end());
159  RCP<const Map> coarseTransferMap = MapFactory::Build(prolongColMap->lib(), INVALID, coarseMapGids(),
160  prolongColMap->getIndexBase(), prolongColMap->getComm());
161 
162  // store map in coarse level
163  if (fineLevel.GetLevelID() == 0) {
164  const std::string mapFactName = pL.get<std::string>("map: factory");
165  RCP<const FactoryBase> mapFact = coarseLevel.GetFactoryManager()->GetFactory(mapFactName);
166  coarseLevel.Set(mapName, coarseTransferMap, mapFact.get());
167  } else
168  coarseLevel.Set(mapName, coarseTransferMap, mapFact_.get());
169 }
170 
171 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
173  const std::string useTheseNspVectors = pL.get<std::string>("nullspace vectors: limit to");
174 
175  // Leave right away, if no limit is prescribed by the user
176  if (useTheseNspVectors == "all" || useTheseNspVectors == "")
177  return -1;
178 
179  // Simplify? Maybe replace by boolean flag "nullspace: exclude rotations"
180  int maxNumProlongCols = -1;
181  if (useTheseNspVectors == "translations")
182  maxNumProlongCols = 1;
183  else
184  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::InvalidArgument, "Unknown subset of nullspace vectors to be used, when performing a map transfer.")
185 
186  return maxNumProlongCols;
187 }
188 
189 } // namespace MueLu
190 
191 #endif /* MUELU_MAPTRANSFERFACTORY_DEF_HPP_ */
ParameterList & setEntry(const std::string &name, U &&entry)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const override
Input.
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
One-liner description of what is happening.
T * get() const
static const NoFactory * get()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< const ParameterList > GetValidParameterList() const override
Input.
int GetLimitOfProlongatorColumns(const ParameterList &pL) const
Get the max number of entries per row of P to be considered for map transfer.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Timer to be used in non-factories.
void Build(Level &fineLevel, Level &coarseLevel) const override
Build an object with this factory.
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:96
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report invalid user entry.
static const RCP< const NoFactory > getRCP()
Static Get() functions.