MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
25 namespace MueLu {
26 
27 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
28 MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MHDRAPFactory()
29  : implicitTranspose_(true) {}
30 
31 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
32 void 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 
57 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
58 void 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 
187 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
188 std::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_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:25
Xpetra_Vector.hpp
MueLu_maxAll
#define MueLu_maxAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:29
GO
GlobalOrdinal GO
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:165
MueLu::MHDRAPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_MHDRAPFactory_def.hpp:58
Xpetra_Matrix.hpp
MueLu_minAll
#define MueLu_minAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:27
MueLu::Statistics2
Print even more statistics.
Definition: MueLu_VerbosityLevel.hpp:35
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
MueLu::MHDRAPFactory::PrintLoadBalancingInfo
static std::string PrintLoadBalancingInfo(const Matrix &Ac, const std::string &msgTag)
Definition: MueLu_MHDRAPFactory_def.hpp:188
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:218
MueLu::MHDRAPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_MHDRAPFactory_def.hpp:32
MueLu_Monitor.hpp
Xpetra_VectorFactory.hpp
MueLu_MHDRAPFactory_decl.hpp