MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_MHDRAPFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MueLu: A package for multigrid based preconditioning
4//
5// Copyright 2012 NTESS and the MueLu contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef MUELU_MHDRAPFACTORY_DEF_HPP
11#define MUELU_MHDRAPFACTORY_DEF_HPP
12
13#include <sstream>
14
15#include <Xpetra_Map.hpp>
16#include <Xpetra_MapFactory.hpp>
17#include <Xpetra_Matrix.hpp>
18#include <Xpetra_CrsMatrixWrap.hpp>
19#include <Xpetra_Vector.hpp>
20#include <Xpetra_VectorFactory.hpp>
21
22#include "MueLu_MHDRAPFactory_decl.hpp"
23#include "MueLu_Monitor.hpp"
24
25namespace MueLu {
26
27template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
30
31template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
32void MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
33 if (implicitTranspose_ == false) {
34 Input(coarseLevel, "R");
35 Input(coarseLevel, "RV");
36 Input(coarseLevel, "RP");
37 Input(coarseLevel, "RM");
38 }
39
40 Input(fineLevel, "A");
41 Input(fineLevel, "A00");
42 Input(fineLevel, "A01");
43 Input(fineLevel, "A02");
44 Input(fineLevel, "A10");
45 Input(fineLevel, "A11");
46 Input(fineLevel, "A12");
47 Input(fineLevel, "A20");
48 Input(fineLevel, "A21");
49 Input(fineLevel, "A22");
50
51 Input(coarseLevel, "P");
52 Input(coarseLevel, "PV");
53 Input(coarseLevel, "PP");
54 Input(coarseLevel, "PM");
55}
56
57template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
58void MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { // FIXME make fineLevel const
59 {
60 FactoryMonitor m(*this, "Computing Ac", coarseLevel);
61
62 //
63 // Inputs: A, P
64 //
65
66 // DEBUG
67 // Teuchos::FancyOStream fout(*GetOStream(Runtime1));
68 // coarseLevel.print(fout,Teuchos::VERB_HIGH);
69
70 RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
71 RCP<Matrix> A00 = Get<RCP<Matrix> >(fineLevel, "A00");
72 RCP<Matrix> A01 = Get<RCP<Matrix> >(fineLevel, "A01");
73 RCP<Matrix> A02 = Get<RCP<Matrix> >(fineLevel, "A02");
74 RCP<Matrix> A10 = Get<RCP<Matrix> >(fineLevel, "A10");
75 RCP<Matrix> A11 = Get<RCP<Matrix> >(fineLevel, "A11");
76 RCP<Matrix> A12 = Get<RCP<Matrix> >(fineLevel, "A12");
77 RCP<Matrix> A20 = Get<RCP<Matrix> >(fineLevel, "A20");
78 RCP<Matrix> A21 = Get<RCP<Matrix> >(fineLevel, "A21");
79 RCP<Matrix> A22 = Get<RCP<Matrix> >(fineLevel, "A22");
80
81 RCP<Matrix> P = Get<RCP<Matrix> >(coarseLevel, "P");
82 RCP<Matrix> PV = Get<RCP<Matrix> >(coarseLevel, "PV");
83 RCP<Matrix> PP = Get<RCP<Matrix> >(coarseLevel, "PP");
84 RCP<Matrix> PM = Get<RCP<Matrix> >(coarseLevel, "PM");
85
86 //
87 // Build Ac = RAP
88 //
89
90 RCP<Matrix> AP;
91 RCP<Matrix> AP00;
92 RCP<Matrix> AP01;
93 RCP<Matrix> AP02;
94 RCP<Matrix> AP10;
95 RCP<Matrix> AP11;
96 RCP<Matrix> AP12;
97 RCP<Matrix> AP20;
98 RCP<Matrix> AP21;
99 RCP<Matrix> AP22;
100
101 {
102 SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
103
104 AP = Utils::Multiply(*A, false, *P, false, AP, GetOStream(Statistics2));
105 AP00 = Utils::Multiply(*A00, false, *PV, false, AP00, GetOStream(Statistics2));
106 AP01 = Utils::Multiply(*A01, false, *PP, false, AP01, GetOStream(Statistics2));
107 AP02 = Utils::Multiply(*A02, false, *PM, false, AP02, GetOStream(Statistics2));
108 AP10 = Utils::Multiply(*A10, false, *PV, false, AP10, GetOStream(Statistics2));
109 AP11 = Utils::Multiply(*A11, false, *PP, false, AP11, GetOStream(Statistics2));
110 AP12 = Utils::Multiply(*A12, false, *PM, false, AP12, GetOStream(Statistics2));
111 AP20 = Utils::Multiply(*A20, false, *PV, false, AP20, GetOStream(Statistics2));
112 AP21 = Utils::Multiply(*A21, false, *PP, false, AP21, GetOStream(Statistics2));
113 AP22 = Utils::Multiply(*A22, false, *PM, false, AP22, GetOStream(Statistics2));
114 }
115
116 RCP<Matrix> Ac;
117 RCP<Matrix> Ac00;
118 RCP<Matrix> Ac01;
119 RCP<Matrix> Ac02;
120 RCP<Matrix> Ac10;
121 RCP<Matrix> Ac11;
122 RCP<Matrix> Ac12;
123 RCP<Matrix> Ac20;
124 RCP<Matrix> Ac21;
125 RCP<Matrix> Ac22;
126
127 if (implicitTranspose_) {
128 SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
129
130 Ac = Utils::Multiply(*P, true, *AP, false, Ac, GetOStream(Statistics2));
131 Ac00 = Utils::Multiply(*PV, true, *AP00, false, Ac00, GetOStream(Statistics2));
132 Ac01 = Utils::Multiply(*PV, true, *AP01, false, Ac01, GetOStream(Statistics2));
133 Ac02 = Utils::Multiply(*PV, true, *AP02, false, Ac02, GetOStream(Statistics2));
134 Ac10 = Utils::Multiply(*PP, true, *AP10, false, Ac10, GetOStream(Statistics2));
135 Ac11 = Utils::Multiply(*PP, true, *AP11, false, Ac11, GetOStream(Statistics2));
136 Ac12 = Utils::Multiply(*PP, true, *AP12, false, Ac12, GetOStream(Statistics2));
137 Ac20 = Utils::Multiply(*PM, true, *AP20, false, Ac20, GetOStream(Statistics2));
138 Ac21 = Utils::Multiply(*PM, true, *AP21, false, Ac21, GetOStream(Statistics2));
139 Ac22 = Utils::Multiply(*PM, true, *AP22, false, Ac22, GetOStream(Statistics2));
140
141 } else {
142 SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
143
144 RCP<Matrix> R = Get<RCP<Matrix> >(coarseLevel, "R");
145 RCP<Matrix> RV = Get<RCP<Matrix> >(coarseLevel, "RV");
146 RCP<Matrix> RP = Get<RCP<Matrix> >(coarseLevel, "RP");
147 RCP<Matrix> RM = Get<RCP<Matrix> >(coarseLevel, "RM");
148
149 Ac = Utils::Multiply(*R, false, *AP, false, Ac, GetOStream(Statistics2));
150 Ac00 = Utils::Multiply(*RV, false, *AP00, false, Ac00, GetOStream(Statistics2));
151 Ac01 = Utils::Multiply(*RV, false, *AP01, false, Ac01, GetOStream(Statistics2));
152 Ac02 = Utils::Multiply(*RV, false, *AP02, false, Ac02, GetOStream(Statistics2));
153 Ac10 = Utils::Multiply(*RP, false, *AP10, false, Ac10, GetOStream(Statistics2));
154 Ac11 = Utils::Multiply(*RP, false, *AP11, false, Ac11, GetOStream(Statistics2));
155 Ac12 = Utils::Multiply(*RP, false, *AP12, false, Ac12, GetOStream(Statistics2));
156 Ac20 = Utils::Multiply(*RM, false, *AP20, false, Ac20, GetOStream(Statistics2));
157 Ac21 = Utils::Multiply(*RM, false, *AP21, false, Ac21, GetOStream(Statistics2));
158 Ac22 = Utils::Multiply(*RM, false, *AP22, false, Ac22, GetOStream(Statistics2));
159 }
160 // FINISHED MAKING COARSE BLOCKS
161
162 Set(coarseLevel, "A", Ac);
163 Set(coarseLevel, "A00", Ac00);
164 Set(coarseLevel, "A01", Ac01);
165 Set(coarseLevel, "A02", Ac02);
166 Set(coarseLevel, "A10", Ac10);
167 Set(coarseLevel, "A11", Ac11);
168 Set(coarseLevel, "A12", Ac12);
169 Set(coarseLevel, "A20", Ac20);
170 Set(coarseLevel, "A21", Ac21);
171 Set(coarseLevel, "A22", Ac22);
172 }
173}
174
175/*
176 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
177 std::string MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PerfUtils::PrintMatrixInfo(const Matrix & Ac, const std::string & msgTag) {
178 std::stringstream ss(std::stringstream::out);
179 ss << msgTag
180 << " # global rows = " << Ac.getGlobalNumRows()
181 << ", estim. global nnz = " << Ac.getGlobalNumEntries()
182 << std::endl;
183 return ss.str();
184 }
185*/
186
187template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
188std::string MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PrintLoadBalancingInfo(const Matrix &Ac, const std::string &msgTag) {
189 std::stringstream ss(std::stringstream::out);
190
191 // TODO: provide a option to skip this (to avoid global communication)
192 // TODO: skip if nproc == 1
193
194 // nonzero imbalance
195 size_t numMyNnz = Ac.getLocalNumEntries();
196 GO maxNnz, minNnz;
197 RCP<const Teuchos::Comm<int> > comm = Ac.getRowMap()->getComm();
198 MueLu_maxAll(comm, (GO)numMyNnz, maxNnz);
199 // min nnz over all proc (disallow any processors with 0 nnz)
200 MueLu_minAll(comm, (GO)((numMyNnz > 0) ? numMyNnz : maxNnz), minNnz);
201 double imbalance = ((double)maxNnz) / minNnz;
202
203 size_t numMyRows = Ac.getLocalNumRows();
204 // Check whether Ac is spread over more than one process.
205 GO numActiveProcesses = 0;
206 MueLu_sumAll(comm, (GO)((numMyRows > 0) ? 1 : 0), numActiveProcesses);
207
208 // min, max, and avg # rows per proc
209 GO minNumRows, maxNumRows;
210 double avgNumRows;
211 MueLu_maxAll(comm, (GO)numMyRows, maxNumRows);
212 MueLu_minAll(comm, (GO)((numMyRows > 0) ? numMyRows : maxNumRows), minNumRows);
213 assert(numActiveProcesses > 0);
214 avgNumRows = Ac.getGlobalNumRows() / numActiveProcesses;
215
216 ss << msgTag << " # processes with rows = " << numActiveProcesses << std::endl;
217 ss << msgTag << " min # rows per proc = " << minNumRows << ", max # rows per proc = " << maxNumRows << ", avg # rows per proc = " << avgNumRows << std::endl;
218 ss << msgTag << " nonzero imbalance = " << imbalance << std::endl;
219
220 return ss.str();
221}
222
223} // namespace MueLu
224
225#define MUELU_MHDRAPFACTORY_SHORT
226#endif // MUELU_MHDRAPFACTORY_DEF_HPP
MueLu_MHDRAPFactory_decl.hpp
MueLu_Monitor.hpp
MueLu_maxAll
#define MueLu_maxAll(rcpComm, in, out)
Definition MueLu_PerfUtils_decl.hpp:29
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition MueLu_PerfUtils_decl.hpp:25
MueLu_minAll
#define MueLu_minAll(rcpComm, in, out)
Definition MueLu_PerfUtils_decl.hpp:27
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition MueLu_Monitor.hpp:165
MueLu::Factory::Input
void Input(Level &level, const std::string &varName) const
Definition MueLu_Factory.cpp:72
MueLu::Factory::Get
T Get(Level &level, const std::string &varName) const
Definition MueLu_Factory.hpp:69
MueLu::Factory::Set
void Set(Level &level, const std::string &varName, const T &data) const
Definition MueLu_Factory.hpp:80
MueLu::Level
Class that holds all level-specific information.
Definition MueLu_Level.hpp:63
MueLu::MHDRAPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition MueLu_MHDRAPFactory_def.hpp:58
MueLu::MHDRAPFactory::PrintLoadBalancingInfo
static std::string PrintLoadBalancingInfo(const Matrix &Ac, const std::string &msgTag)
Definition MueLu_MHDRAPFactory_def.hpp:188
MueLu::MHDRAPFactory::implicitTranspose_
bool implicitTranspose_
If true, the action of the restriction operator action is implicitly defined by the transpose of the ...
Definition MueLu_MHDRAPFactory_decl.hpp:76
MueLu::MHDRAPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition MueLu_MHDRAPFactory_def.hpp:32
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition MueLu_Monitor.hpp:218
MueLu::VerboseObject::GetOStream
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Definition MueLu_VerboseObject.cpp:72
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:42
MueLu::Statistics2
@ Statistics2
Print even more statistics.
Definition MueLu_VerbosityLevel.hpp:35