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 // ***********************************************************************
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_MHDRAPFACTORY_DEF_HPP
47 #define MUELU_MHDRAPFACTORY_DEF_HPP
48 
49 #include <sstream>
50 
51 #include <Xpetra_Map.hpp>
52 #include <Xpetra_MapFactory.hpp>
53 #include <Xpetra_Matrix.hpp>
54 #include <Xpetra_CrsMatrixWrap.hpp>
55 #include <Xpetra_Vector.hpp>
56 #include <Xpetra_VectorFactory.hpp>
57 
58 #include "MueLu_MHDRAPFactory_decl.hpp"
59 #include "MueLu_Utilities.hpp"
60 #include "MueLu_Monitor.hpp"
61 
62 namespace MueLu {
63 
64  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MHDRAPFactory()
66  : implicitTranspose_(true) { }
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  void MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
70 
71  if (implicitTranspose_ == false)
72  {
73  Input(coarseLevel, "R" );
74  Input(coarseLevel, "RV");
75  Input(coarseLevel, "RP");
76  Input(coarseLevel, "RM");
77  }
78 
79  Input(fineLevel, "A" );
80  Input(fineLevel, "A00");
81  Input(fineLevel, "A01");
82  Input(fineLevel, "A02");
83  Input(fineLevel, "A10");
84  Input(fineLevel, "A11");
85  Input(fineLevel, "A12");
86  Input(fineLevel, "A20");
87  Input(fineLevel, "A21");
88  Input(fineLevel, "A22");
89 
90  Input(coarseLevel, "P" );
91  Input(coarseLevel, "PV");
92  Input(coarseLevel, "PP");
93  Input(coarseLevel, "PM");
94 
95  }
96 
97  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98  void MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { // FIXME make fineLevel const
99  {
100  FactoryMonitor m(*this, "Computing Ac", coarseLevel);
101 
102  //
103  // Inputs: A, P
104  //
105 
106  //DEBUG
107  //Teuchos::FancyOStream fout(*GetOStream(Runtime1));
108  //coarseLevel.print(fout,Teuchos::VERB_HIGH);
109 
110 
111  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A" );
112  RCP<Matrix> A00 = Get< RCP<Matrix> >(fineLevel, "A00");
113  RCP<Matrix> A01 = Get< RCP<Matrix> >(fineLevel, "A01");
114  RCP<Matrix> A02 = Get< RCP<Matrix> >(fineLevel, "A02");
115  RCP<Matrix> A10 = Get< RCP<Matrix> >(fineLevel, "A10");
116  RCP<Matrix> A11 = Get< RCP<Matrix> >(fineLevel, "A11");
117  RCP<Matrix> A12 = Get< RCP<Matrix> >(fineLevel, "A12");
118  RCP<Matrix> A20 = Get< RCP<Matrix> >(fineLevel, "A20");
119  RCP<Matrix> A21 = Get< RCP<Matrix> >(fineLevel, "A21");
120  RCP<Matrix> A22 = Get< RCP<Matrix> >(fineLevel, "A22");
121 
122  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P" );
123  RCP<Matrix> PV = Get< RCP<Matrix> >(coarseLevel, "PV");
124  RCP<Matrix> PP = Get< RCP<Matrix> >(coarseLevel, "PP");
125  RCP<Matrix> PM = Get< RCP<Matrix> >(coarseLevel, "PM");
126 
127  //
128  // Build Ac = RAP
129  //
130 
131  RCP<Matrix> AP;
132  RCP<Matrix> AP00;
133  RCP<Matrix> AP01;
134  RCP<Matrix> AP02;
135  RCP<Matrix> AP10;
136  RCP<Matrix> AP11;
137  RCP<Matrix> AP12;
138  RCP<Matrix> AP20;
139  RCP<Matrix> AP21;
140  RCP<Matrix> AP22;
141 
142  {
143  SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
144 
145  AP = Utils::Multiply(*A , false, *P , false, AP, GetOStream(Statistics2));
146  AP00 = Utils::Multiply(*A00, false, *PV, false, AP00, GetOStream(Statistics2));
147  AP01 = Utils::Multiply(*A01, false, *PP, false, AP01, GetOStream(Statistics2));
148  AP02 = Utils::Multiply(*A02, false, *PM, false, AP02, GetOStream(Statistics2));
149  AP10 = Utils::Multiply(*A10, false, *PV, false, AP10, GetOStream(Statistics2));
150  AP11 = Utils::Multiply(*A11, false, *PP, false, AP11, GetOStream(Statistics2));
151  AP12 = Utils::Multiply(*A12, false, *PM, false, AP12, GetOStream(Statistics2));
152  AP20 = Utils::Multiply(*A20, false, *PV, false, AP20, GetOStream(Statistics2));
153  AP21 = Utils::Multiply(*A21, false, *PP, false, AP21, GetOStream(Statistics2));
154  AP22 = Utils::Multiply(*A22, false, *PM, false, AP22, GetOStream(Statistics2));
155  }
156 
157  RCP<Matrix> Ac;
158  RCP<Matrix> Ac00;
159  RCP<Matrix> Ac01;
160  RCP<Matrix> Ac02;
161  RCP<Matrix> Ac10;
162  RCP<Matrix> Ac11;
163  RCP<Matrix> Ac12;
164  RCP<Matrix> Ac20;
165  RCP<Matrix> Ac21;
166  RCP<Matrix> Ac22;
167 
168  if (implicitTranspose_)
169  {
170  SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
171 
172  Ac = Utils::Multiply(*P , true, *AP , false, Ac, GetOStream(Statistics2));
173  Ac00 = Utils::Multiply(*PV, true, *AP00, false, Ac00, GetOStream(Statistics2));
174  Ac01 = Utils::Multiply(*PV, true, *AP01, false, Ac01, GetOStream(Statistics2));
175  Ac02 = Utils::Multiply(*PV, true, *AP02, false, Ac02, GetOStream(Statistics2));
176  Ac10 = Utils::Multiply(*PP, true, *AP10, false, Ac10, GetOStream(Statistics2));
177  Ac11 = Utils::Multiply(*PP, true, *AP11, false, Ac11, GetOStream(Statistics2));
178  Ac12 = Utils::Multiply(*PP, true, *AP12, false, Ac12, GetOStream(Statistics2));
179  Ac20 = Utils::Multiply(*PM, true, *AP20, false, Ac20, GetOStream(Statistics2));
180  Ac21 = Utils::Multiply(*PM, true, *AP21, false, Ac21, GetOStream(Statistics2));
181  Ac22 = Utils::Multiply(*PM, true, *AP22, false, Ac22, GetOStream(Statistics2));
182 
183  }
184  else
185  {
186 
187  SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
188 
189  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R" );
190  RCP<Matrix> RV = Get< RCP<Matrix> >(coarseLevel, "RV");
191  RCP<Matrix> RP = Get< RCP<Matrix> >(coarseLevel, "RP");
192  RCP<Matrix> RM = Get< RCP<Matrix> >(coarseLevel, "RM");
193 
194  Ac = Utils::Multiply(*R , false, *AP , false, Ac, GetOStream(Statistics2));
195  Ac00 = Utils::Multiply(*RV, false, *AP00, false, Ac00, GetOStream(Statistics2));
196  Ac01 = Utils::Multiply(*RV, false, *AP01, false, Ac01, GetOStream(Statistics2));
197  Ac02 = Utils::Multiply(*RV, false, *AP02, false, Ac02, GetOStream(Statistics2));
198  Ac10 = Utils::Multiply(*RP, false, *AP10, false, Ac10, GetOStream(Statistics2));
199  Ac11 = Utils::Multiply(*RP, false, *AP11, false, Ac11, GetOStream(Statistics2));
200  Ac12 = Utils::Multiply(*RP, false, *AP12, false, Ac12, GetOStream(Statistics2));
201  Ac20 = Utils::Multiply(*RM, false, *AP20, false, Ac20, GetOStream(Statistics2));
202  Ac21 = Utils::Multiply(*RM, false, *AP21, false, Ac21, GetOStream(Statistics2));
203  Ac22 = Utils::Multiply(*RM, false, *AP22, false, Ac22, GetOStream(Statistics2));
204 
205  }
206  // FINISHED MAKING COARSE BLOCKS
207 
208  Set(coarseLevel, "A" , Ac );
209  Set(coarseLevel, "A00", Ac00);
210  Set(coarseLevel, "A01", Ac01);
211  Set(coarseLevel, "A02", Ac02);
212  Set(coarseLevel, "A10", Ac10);
213  Set(coarseLevel, "A11", Ac11);
214  Set(coarseLevel, "A12", Ac12);
215  Set(coarseLevel, "A20", Ac20);
216  Set(coarseLevel, "A21", Ac21);
217  Set(coarseLevel, "A22", Ac22);
218 
219 
220  }
221 
222 
223  }
224 
225 
226 /*
227  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
228  std::string MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PerfUtils::PrintMatrixInfo(const Matrix & Ac, const std::string & msgTag) {
229  std::stringstream ss(std::stringstream::out);
230  ss << msgTag
231  << " # global rows = " << Ac.getGlobalNumRows()
232  << ", estim. global nnz = " << Ac.getGlobalNumEntries()
233  << std::endl;
234  return ss.str();
235  }
236 */
237 
238  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
239  std::string MHDRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PrintLoadBalancingInfo(const Matrix & Ac, const std::string & msgTag) {
240  std::stringstream ss(std::stringstream::out);
241 
242  // TODO: provide a option to skip this (to avoid global communication)
243  // TODO: skip if nproc == 1
244 
245  //nonzero imbalance
246  size_t numMyNnz = Ac.getNodeNumEntries();
247  GO maxNnz, minNnz;
248  RCP<const Teuchos::Comm<int> > comm = Ac.getRowMap()->getComm();
249  MueLu_maxAll(comm,(GO)numMyNnz,maxNnz);
250  //min nnz over all proc (disallow any processors with 0 nnz)
251  MueLu_minAll(comm, (GO)((numMyNnz > 0) ? numMyNnz : maxNnz), minNnz);
252  double imbalance = ((double) maxNnz) / minNnz;
253 
254  size_t numMyRows = Ac.getNodeNumRows();
255  //Check whether Ac is spread over more than one process.
256  GO numActiveProcesses=0;
257  MueLu_sumAll(comm, (GO)((numMyRows > 0) ? 1 : 0), numActiveProcesses);
258 
259  //min, max, and avg # rows per proc
260  GO minNumRows, maxNumRows;
261  double avgNumRows;
262  MueLu_maxAll(comm, (GO)numMyRows, maxNumRows);
263  MueLu_minAll(comm, (GO)((numMyRows > 0) ? numMyRows : maxNumRows), minNumRows);
264  assert(numActiveProcesses > 0);
265  avgNumRows = Ac.getGlobalNumRows() / numActiveProcesses;
266 
267  ss << msgTag << " # processes with rows = " << numActiveProcesses << std::endl;
268  ss << msgTag << " min # rows per proc = " << minNumRows << ", max # rows per proc = " << maxNumRows << ", avg # rows per proc = " << avgNumRows << std::endl;
269  ss << msgTag << " nonzero imbalance = " << imbalance << std::endl;
270 
271  return ss.str();
272  }
273 
274 
275 } //namespace MueLu
276 
277 #define MUELU_MHDRAPFACTORY_SHORT
278 #endif // MUELU_MHDRAPFACTORY_DEF_HPP
279 
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:61
Xpetra_Vector.hpp
MueLu_maxAll
#define MueLu_maxAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:65
GO
GlobalOrdinal GO
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
MueLu::MHDRAPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_MHDRAPFactory_def.hpp:98
Xpetra_Matrix.hpp
MueLu_minAll
#define MueLu_minAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:63
MueLu::Statistics2
Print even more statistics.
Definition: MueLu_VerbosityLevel.hpp:72
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::MHDRAPFactory::PrintLoadBalancingInfo
static std::string PrintLoadBalancingInfo(const Matrix &Ac, const std::string &msgTag)
Definition: MueLu_MHDRAPFactory_def.hpp:239
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:261
Xpetra_MapFactory.hpp
MueLu::MHDRAPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_MHDRAPFactory_def.hpp:69
MueLu_Monitor.hpp
Xpetra_VectorFactory.hpp
Xpetra_Map.hpp
Xpetra_CrsMatrixWrap.hpp
MueLu_MHDRAPFactory_decl.hpp