MueLu Version of the Day
Loading...
Searching...
No Matches
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
59#include "MueLu_Utilities.hpp"
60#include "MueLu_Monitor.hpp"
61
62namespace MueLu {
63
64 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
67
68 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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.getLocalNumEntries();
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.getLocalNumRows();
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
#define MueLu_maxAll(rcpComm, in, out)
#define MueLu_sumAll(rcpComm, in, out)
#define MueLu_minAll(rcpComm, in, out)
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
static std::string PrintLoadBalancingInfo(const Matrix &Ac, const std::string &msgTag)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.