MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_NullspaceFactory_kokkos_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#ifndef MUELU_NULLSPACEFACTORY_KOKKOS_DEF_HPP
47#define MUELU_NULLSPACEFACTORY_KOKKOS_DEF_HPP
48
50
51#include <Xpetra_Matrix.hpp>
52#include <Xpetra_MultiVectorFactory.hpp>
53
54#include "MueLu_Level.hpp"
55#include "MueLu_Monitor.hpp"
56
57namespace MueLu {
58
59 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
61 RCP<ParameterList> validParamList = rcp(new ParameterList());
62
63 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the fine level matrix (only needed if default null space is generated)");
64 validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the fine level null space");
65
66 // TODO not very elegant.
67 validParamList->set< std::string >("Fine level nullspace", "Nullspace", "Variable name which is used to store null space multi vector on the finest level (default=\"Nullspace\"). For block matrices also \"Nullspace1\" to \"Nullspace9\" are accepted to describe the null space vectors for the (i,i) block (i=1..9).");
68
69 // 1/20/2016: we could add a sublist (e.g. "Nullspaces" which is excluded from parameter validation)
70 validParamList->set< RCP<const FactoryBase> >("Nullspace1", Teuchos::null, "Generating factory of the fine level null space associated with the (1,1) block in your n x n block matrix.");
71 validParamList->set< RCP<const FactoryBase> >("Nullspace2", Teuchos::null, "Generating factory of the fine level null space associated with the (2,2) block in your n x n block matrix.");
72 validParamList->set< RCP<const FactoryBase> >("Nullspace3", Teuchos::null, "Generating factory of the fine level null space associated with the (3,3) block in your n x n block matrix.");
73 validParamList->set< RCP<const FactoryBase> >("Nullspace4", Teuchos::null, "Generating factory of the fine level null space associated with the (4,4) block in your n x n block matrix.");
74 validParamList->set< RCP<const FactoryBase> >("Nullspace5", Teuchos::null, "Generating factory of the fine level null space associated with the (5,5) block in your n x n block matrix.");
75 validParamList->set< RCP<const FactoryBase> >("Nullspace6", Teuchos::null, "Generating factory of the fine level null space associated with the (6,6) block in your n x n block matrix.");
76 validParamList->set< RCP<const FactoryBase> >("Nullspace7", Teuchos::null, "Generating factory of the fine level null space associated with the (7,7) block in your n x n block matrix.");
77 validParamList->set< RCP<const FactoryBase> >("Nullspace8", Teuchos::null, "Generating factory of the fine level null space associated with the (8,8) block in your n x n block matrix.");
78 validParamList->set< RCP<const FactoryBase> >("Nullspace9", Teuchos::null, "Generating factory of the fine level null space associated with the (9,9) block in your n x n block matrix.");
79
80 return validParamList;
81 }
82
83 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
85 const ParameterList& pL = GetParameterList();
86 std::string nspName = pL.get<std::string>("Fine level nullspace");
87
88 // only request "A" in DeclareInput if
89 // 1) there is not nspName (e.g. "Nullspace") is available in Level, AND
90 // 2) it is the finest level (i.e. LevelID == 0)
91 if (currentLevel.IsAvailable(nspName, NoFactory::get()) == false && currentLevel.GetLevelID() == 0)
92 Input(currentLevel, "A");
93
94 if (currentLevel.GetLevelID() != 0) {
95 // validate nullspaceFact_
96 // 1) The factory for "Nullspace" (or nspName) must not be Teuchos::null, since the default factory
97 // for "Nullspace" is a NullspaceFactory
98 // 2) The factory for "Nullspace" (or nspName) must be a TentativePFactory or any other TwoLevelFactoryBase derived object
99 // which generates the variable "Nullspace" as output
100 TEUCHOS_TEST_FOR_EXCEPTION(GetFactory(nspName).is_null(), Exceptions::RuntimeError,
101 "MueLu::NullspaceFactory::DeclareInput(): You must declare an existing factory which "
102 "produces the variable \"Nullspace\" in the NullspaceFactory (e.g. a TentativePFactory).");
103 currentLevel.DeclareInput("Nullspace", GetFactory(nspName).get(), this); /* ! "Nullspace" and nspName mismatch possible here */
104 }
105 }
106
107 template<class NullspaceType, class LO>
109 private:
110 NullspaceType nullspace;
112 typedef typename NullspaceType::value_type SC;
113 typedef Kokkos::ArithTraits<SC> ATS;
114
115 public:
116 NullspaceFunctor(NullspaceType nullspace_, LO numPDEs_) :
117 nullspace(nullspace_),
118 numPDEs(numPDEs_)
119 { }
120
121 KOKKOS_INLINE_FUNCTION
122 void operator()(const LO j) const {
123 SC one = ATS::one();
124 for (LO i = 0; i < numPDEs; i++)
125 nullspace(j*numPDEs + i, i) = one;
126 }
127 };
128 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
130 FactoryMonitor m(*this, "Nullspace factory", currentLevel);
131
132 RCP<MultiVector> nullspace;
133
134 //TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.GetLevelID() != 0, Exceptions::RuntimeError, "MueLu::NullspaceFactory::Build(): NullspaceFactory can be used for finest level (LevelID == 0) only.");
135 const ParameterList& pL = GetParameterList();
136 std::string nspName = pL.get<std::string>("Fine level nullspace");
137
138 if (currentLevel.GetLevelID() == 0) {
139 if (currentLevel.IsAvailable(nspName, NoFactory::get())) {
140 // When a fine nullspace have already been defined by user using Set("Nullspace", ...) or
141 // Set("Nullspace1", ...), we use it.
142 nullspace = currentLevel.Get< RCP<MultiVector> >(nspName, NoFactory::get());
143 GetOStream(Runtime1) << "Use user-given nullspace " << nspName << ":"
144 << " nullspace dimension=" << nullspace->getNumVectors()
145 << " nullspace length=" << nullspace->getGlobalLength() << std::endl;
146
147 } else {
148 // "Nullspace" (nspName) is not available
149 auto A = Get< RCP<Matrix> >(currentLevel, "A");
150
151 // Determine numPDEs
152 LO numPDEs = 1;
153 if (A->IsView("stridedMaps") == true) {
154 Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
155 TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<const StridedMap>(A->getRowMap()).is_null(), Exceptions::BadCast,
156 "MueLu::CoalesceFactory::Build: cast to strided row map failed.");
157 numPDEs = rcp_dynamic_cast<const StridedMap>(A->getRowMap())->getFixedBlockSize();
158 oldView = A->SwitchToView(oldView);
159 }
160
161 GetOStream(Runtime1) << "Generating canonical nullspace: dimension = " << numPDEs << std::endl;
162
163 nullspace = MultiVectorFactory::Build(A->getDomainMap(), numPDEs);
164 auto nullspaceView = nullspace->getDeviceLocalView(Xpetra::Access::OverwriteAll);
165
166 int numBlocks = nullspace->getLocalLength() / numPDEs;
167
168 NullspaceFunctor<decltype(nullspaceView), LO> nullspaceFunctor(nullspaceView, numPDEs);
169 Kokkos::parallel_for("MueLu:NullspaceF:Build:for", range_type(0,numBlocks), nullspaceFunctor);
170 // TODO extend null space factory for blocked multi vectors
171 }
172
173 } else {
174 // On coarser levels always use "Nullspace" as variable name, since it is expected by
175 // tentative P factory to be "Nullspace"
176 nullspace = currentLevel.Get< RCP<MultiVector> >("Nullspace", GetFactory(nspName).get()); // NOTE: "Nullspace" and nspName mismatch possible here
177 }
178
179 // provide "Nullspace" variable on current level (used by TentativePFactory)
180 Set(currentLevel, "Nullspace", nullspace);
181
182 } // Build
183
184} //namespace MueLu
185
186#endif // MUELU_NULLSPACEFACTORY_KOKKOS_DEF_HPP
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.
Definition: MueLu_Level.hpp:99
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()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
NullspaceFunctor(NullspaceType nullspace_, LO numPDEs_)
KOKKOS_INLINE_FUNCTION void operator()(const LO j) const
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose)