MueLu Version of the Day
MueLu_PerfUtils_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 NodeT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DIScalarLAIMED. IN Node EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NodeT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GlobalOrdinalODS 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_PERFUTILS_DEF_HPP
47#define MUELU_PERFUTILS_DEF_HPP
48
49#include <algorithm>
50#include <string>
51
52#ifdef HAVE_MPI
53#include <Teuchos_CommHelpers.hpp>
54#endif
55
56#include <Xpetra_Export.hpp>
57#include <Xpetra_Import.hpp>
58#include <Xpetra_Matrix.hpp>
59#include <Xpetra_CrsMatrixWrap.hpp>
60
61#include "MueLu_PerfUtils_decl.hpp"
62
63//#include "MueLu_Utilities.hpp"
64
65namespace MueLu {
66
67 template<class Type>
68 void calculateStats(Type& minVal, Type& maxVal, double& avgVal, double& devVal, int& minProc, int& maxProc, const RCP<const Teuchos::Comm<int> >& comm, int numActiveProcs, const Type& v) {
69
70 Type sumVal, sum2Val, v2 = v*v;
71 double zero = Teuchos::ScalarTraits<double>::zero();
72
73 MueLu_sumAll(comm, v, sumVal);
74 MueLu_sumAll(comm, v2, sum2Val);
75 MueLu_minAll(comm, v, minVal);
76 MueLu_maxAll(comm, v, maxVal);
77
78 int w;
79 w = (minVal == v) ? comm->getRank() : -1;
80 MueLu_maxAll(comm, w, maxProc);
81 w = (maxVal == v) ? comm->getRank() : -1;
82 MueLu_maxAll(comm, w, minProc);
83
84 avgVal = (numActiveProcs > 0 ? (as<double>(Teuchos::ScalarTraits<Type>::real(sumVal)) / numActiveProcs) : zero);
85 devVal = (numActiveProcs > 1 ? sqrt((as<double>(Teuchos::ScalarTraits<Type>::real(sum2Val - sumVal*avgVal)))/(numActiveProcs-1)) : zero);
86 }
87
88 template<class Type>
89 std::string stringStats(const RCP<const Teuchos::Comm<int> >& comm, int numActiveProcs, const Type& v, RCP<ParameterList> paramList = Teuchos::null) {
90 Type minVal, maxVal;
91 double avgVal, devVal;
92 int minProc, maxProc;
93 calculateStats<Type>(minVal, maxVal, avgVal, devVal, minProc, maxProc, comm, numActiveProcs, v);
94
95 const double zero = Teuchos::ScalarTraits<double>::zero();
96 const double one = Teuchos::ScalarTraits<double>::one();
97 std::ostringstream buf;
98 buf << std::fixed;
99 if ((avgVal != zero) && (paramList.is_null() || !paramList->isParameter("print abs") || paramList->get<bool>("print abs") == false)) {
100 double relDev = (devVal/avgVal)*100;
101 double relMin = (as<double>(Teuchos::ScalarTraits<Type>::real(minVal))/avgVal-one)*100;
102 double relMax = (as<double>(Teuchos::ScalarTraits<Type>::real(maxVal))/avgVal-one)*100;
103 buf << "avg = " << std::scientific << std::setw(10) << std::setprecision(2) << avgVal << ", "
104 << "dev = " << std::fixed << std::setw(6) << std::setprecision(1) << relDev << "%, "
105 << "min = " << std::fixed << std::setw(7) << std::setprecision(1) << std::setw(7) << relMin << "%"
106 << " (" << std::scientific << std::setw(10) << std::setprecision(2) << minVal << " on " << std::fixed << std::setw(4) << minProc << "), "
107 << "max = " << std::fixed << std::setw(7) << std::setprecision(1) << relMax << "%"
108 << " (" << std::scientific << std::setw(10) << std::setprecision(2) << maxVal << " on " << std::fixed << std::setw(4) << maxProc << ")";
109 } else {
110 double relDev = (avgVal != zero ? (devVal/avgVal)*100 : zero);
111 buf << "avg = " << std::scientific << std::setw(10) << std::setprecision(2) << avgVal << ", "
112 << "dev = " << std::fixed << std::setw(6) << std::setprecision(1) << relDev << "%, "
113 << "min = " << std::scientific << std::setw(10) << std::setprecision(2) << minVal
114 << " (on " << std::fixed << std::setw(4) << minProc << "), "
115 << "max = " << std::scientific << std::setw(10) << std::setprecision(2) << maxVal
116 << " (on " << std::fixed << std::setw(4) << maxProc << ")";
117 }
118 return buf.str();
119 }
120
121 template<class Map>
122 bool cmp_less(typename Map::value_type& v1, typename Map::value_type& v2) {
123 return v1.second < v2.second;
124 }
125
126 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
127 std::string PerfUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PrintMatrixInfo(const Matrix& A, const std::string& msgTag, RCP<const ParameterList> params) {
128 if (!CheckMatrix(A))
129 return "";
130
131 typedef Xpetra::global_size_t global_size_t;
132
133 std::ostringstream ss;
134
135 ss << msgTag << " size = " << A.getGlobalNumRows() << " x " << A.getGlobalNumCols();
136 if(A.haveGlobalConstants())
137 ss << ", nnz = " << A.getGlobalNumEntries();
138 ss << std::endl;
139
140 if (params.is_null())
141 return ss.str();
142
143 bool printLoadBalanceInfo = false, printCommInfo = false, printEntryStats = false;
144 if (params->isParameter("printLoadBalancingInfo") && params->get<bool>("printLoadBalancingInfo"))
145 printLoadBalanceInfo = true;
146 if (params->isParameter("printCommInfo") && params->get<bool>("printCommInfo"))
147 printCommInfo = true;
148 if (params->isParameter("printEntryStats") && params->get<bool>("printEntryStats"))
149 printEntryStats = true;
150
151 if (!printLoadBalanceInfo && !printCommInfo && !printEntryStats)
152 return ss.str();
153
154 RCP<const Import> importer = A.getCrsGraph()->getImporter();
155 RCP<const Export> exporter = A.getCrsGraph()->getExporter();
156
157 size_t numMyNnz = A.getNodeNumEntries(), numMyRows = A.getNodeNumRows();
158
159 // Create communicator only for active processes
160 RCP<const Teuchos::Comm<int> > origComm = A.getRowMap()->getComm();
161 bool activeProc = true;
162 int numProc = origComm->getSize();
163 int numActiveProcs = 0;
164#ifdef HAVE_MPI
165 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(origComm);
166 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
167
168 std::vector<size_t> numRowsPerProc(numProc);
169 Teuchos::gatherAll(*origComm, 1, &numMyRows, numProc, &numRowsPerProc[0]);
170
171 int root = 0;
172 bool rootFlag = true;
173 for (int i = 0; i < numProc; i++) {
174 if (numRowsPerProc[i]) {
175 ++numActiveProcs;
176 if(rootFlag) {
177 root = i;
178 rootFlag = false;
179 }
180 }
181 }
182
183 if(numMyRows == 0) {activeProc = false; numMyNnz = 0;} // Reset numMyNnz to avoid adding it up in reduceAll
184#else
185 if(numMyRows == 0) {
186 //FIXME JJH 10-May-2017 Is there any case in serial where numMyRows would be zero?
187 // Reset numMyNnz to avoid adding it up in reduceAll
188 numActiveProcs = 0;
189 activeProc = false;
190 numMyNnz = 0;
191 } else {
192 numActiveProcs = 1;
193 }
194#endif
195
196 std::string outstr;
197 ParameterList absList;
198 absList.set("print abs", true);
199
200 RCP<const Matrix> rcpA = rcpFromRef(A);
201 RCP<const CrsMatrixWrap> crsWrapA = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(rcpA);
202 RCP<const CrsMatrix> crsA;
203 if (!crsWrapA.is_null())
204 crsA = crsWrapA->getCrsMatrix();
205 if (printEntryStats && !crsA.is_null()) {
206 typedef Teuchos::ScalarTraits<Scalar> STS;
207 typedef typename STS::magnitudeType magnitudeType;
208 typedef Teuchos::ScalarTraits<magnitudeType> MTS;
209 ArrayRCP<const size_t> rowptr_RCP;
210 ArrayRCP<const LocalOrdinal> colind_RCP;
211 ArrayRCP<const Scalar> vals_RCP;
212 ArrayRCP<size_t> offsets_RCP;
213 ArrayView<const size_t> rowptr;
214 ArrayView<const Scalar> vals;
215 ArrayView<size_t> offsets;
216
217 crsA->getAllValues(rowptr_RCP, colind_RCP, vals_RCP);
218 crsA->getLocalDiagOffsets(offsets_RCP);
219 rowptr = rowptr_RCP();
220 vals = vals_RCP();
221 offsets = offsets_RCP();
222
223 Scalar val, minVal, maxVal;
224 magnitudeType absVal, minAbsVal, maxAbsVal;
225 {
226 minVal = STS::rmax();
227 maxVal = STS::rmin();
228 minAbsVal = MTS::rmax();
229 maxAbsVal = MTS::zero();
230
231 for (int i = 0; i < offsets.size(); i++) {
232 val = vals[rowptr[i]+offsets[i]];
233 if (STS::real(val) < STS::real(minVal))
234 minVal = val;
235 if (STS::real(val) > STS::real(maxVal))
236 maxVal = val;
237 absVal = STS::magnitude(val);
238 minAbsVal = std::min(minAbsVal, absVal);
239 maxAbsVal = std::max(maxAbsVal, absVal);
240 }
241
242 ss << msgTag << " diag min : " << stringStats<Scalar>(origComm, numActiveProcs, minVal) << std::endl;
243 ss << msgTag << " diag max : " << stringStats<Scalar>(origComm, numActiveProcs, maxVal) << std::endl;
244 ss << msgTag << " abs(diag) min : " << stringStats<Scalar>(origComm, numActiveProcs, minAbsVal) << std::endl;
245 ss << msgTag << " abs(diag) max : " << stringStats<Scalar>(origComm, numActiveProcs, maxAbsVal) << std::endl;
246 }
247
248 {
249 minVal = STS::rmax();
250 maxVal = STS::rmin();
251 minAbsVal = MTS::rmax();
252 maxAbsVal = MTS::zero();
253
254 for (int i = 0; i < vals.size(); i++) {
255 val = vals[i];
256 if (STS::real(val) < STS::real(minVal))
257 minVal = val;
258 if (STS::real(val) > STS::real(maxVal))
259 maxVal = val;
260 absVal = STS::magnitude(val);
261 minAbsVal = std::min(minAbsVal, absVal);
262 maxAbsVal = std::max(maxAbsVal, absVal);
263 }
264
265 ss << msgTag << " entry min : " << stringStats<Scalar>(origComm, numActiveProcs, minVal) << std::endl;
266 ss << msgTag << " entry max : " << stringStats<Scalar>(origComm, numActiveProcs, maxVal) << std::endl;
267 ss << msgTag << " abs(entry) min : " << stringStats<Scalar>(origComm, numActiveProcs, minAbsVal) << std::endl;
268 ss << msgTag << " abs(entry) max : " << stringStats<Scalar>(origComm, numActiveProcs, maxAbsVal) << std::endl;
269 }
270 }
271
272
273 if (printLoadBalanceInfo) {
274 ss << msgTag << " Load balancing info" << std::endl;
275 ss << msgTag << " # active processes: " << numActiveProcs << "/" << numProc << std::endl;
276 ss << msgTag << " # rows per proc : " << stringStats<global_size_t>(origComm, numActiveProcs, numMyRows) << std::endl;
277 ss << msgTag << " # nnz per proc : " << stringStats<global_size_t>(origComm, numActiveProcs, numMyNnz) << std::endl;
278 }
279
280 if (printCommInfo && numActiveProcs != 1) {
281 typedef std::map<int,size_t> map_type;
282 map_type neighMap;
283 if (!importer.is_null()) {
284 ArrayView<const int> exportPIDs = importer->getExportPIDs();
285 if (exportPIDs.size())
286 for (int i = 0; i < exportPIDs.size(); i++)
287 neighMap[exportPIDs[i]]++;
288 }
289
290 // Communication volume
291 size_t numExportSend = 0;
292 size_t numImportSend = 0;
293 size_t numMsgs = 0;
294 size_t minMsg = 0;
295 size_t maxMsg = 0;
296
297 if(activeProc) {
298 numExportSend = (!exporter.is_null() ? exporter->getNumExportIDs() : 0);
299 numImportSend = (!importer.is_null() ? importer->getNumExportIDs() : 0);
300 numMsgs = neighMap.size();
301 map_type::const_iterator it = std::min_element(neighMap.begin(), neighMap.end(), cmp_less<map_type>);
302 minMsg = (it != neighMap.end() ? it->second : 0);
303 it = std::max_element(neighMap.begin(), neighMap.end(), cmp_less<map_type>);
304 maxMsg = (it != neighMap.end() ? it->second : 0);
305 }
306
307 ss << msgTag << " Communication info" << std::endl;
308 ss << msgTag << " # num export send : " << stringStats<global_size_t>(origComm, numActiveProcs, numExportSend) << std::endl;
309 ss << msgTag << " # num import send : " << stringStats<global_size_t>(origComm, numActiveProcs, numImportSend) << std::endl;
310 ss << msgTag << " # num msgs : " << stringStats<global_size_t>(origComm, numActiveProcs, numMsgs, rcpFromRef(absList)) << std::endl;
311 ss << msgTag << " # min msg size : " << stringStats<global_size_t>(origComm, numActiveProcs, minMsg) << std::endl;
312 ss << msgTag << " # max msg size : " << stringStats<global_size_t>(origComm, numActiveProcs, maxMsg) << std::endl;
313 }
314
315 outstr = ss.str();
316
317#ifdef HAVE_MPI
318 int strLength = outstr.size();
319 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
320 if (origComm->getRank() != root)
321 outstr.resize(strLength);
322 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
323#endif
324
325 return outstr;
326 }
327
328 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329 std::string PerfUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PrintImporterInfo(RCP<const Import> importer, const std::string& msgTag) {
330
331 typedef Xpetra::global_size_t global_size_t;
332
333 std::ostringstream ss;
334
335 // Create communicator only for active processes
336 RCP<const Teuchos::Comm<int> > origComm = importer->getSourceMap()->getComm();
337 bool activeProc = true;
338 int numActiveProcs = origComm->getSize();
339#ifdef HAVE_MPI
340 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(origComm);
341 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
342 int root = 0;
343#endif
344
345 std::string outstr;
346 ParameterList absList;
347 absList.set("print abs", true);
348
349 typedef std::map<int,size_t> map_type;
350 map_type neighMap;
351 ArrayView<const int> exportPIDs = importer->getExportPIDs();
352 if (exportPIDs.size())
353 for (int i = 0; i < exportPIDs.size(); i++)
354 neighMap[exportPIDs[i]]++;
355
356 // Communication volume
357 size_t numImportSend = 0;
358 size_t numMsgs = 0;
359 size_t minMsg = 0;
360 size_t maxMsg = 0;
361
362 if(activeProc) {
363 numImportSend = importer->getNumExportIDs();
364 numMsgs = neighMap.size();
365 map_type::const_iterator it = std::min_element(neighMap.begin(), neighMap.end(), cmp_less<map_type>);
366 minMsg = (it != neighMap.end() ? it->second : 0);
367 it = std::max_element(neighMap.begin(), neighMap.end(), cmp_less<map_type>);
368 maxMsg = (it != neighMap.end() ? it->second : 0);
369 }
370
371 ss << msgTag << " Communication info" << std::endl;
372 ss << msgTag << " # num import send : " << stringStats<global_size_t>(origComm, numActiveProcs, numImportSend) << std::endl;
373 ss << msgTag << " # num msgs : " << stringStats<global_size_t>(origComm, numActiveProcs, numMsgs, rcpFromRef(absList)) << std::endl;
374 ss << msgTag << " # min msg size : " << stringStats<global_size_t>(origComm, numActiveProcs, minMsg) << std::endl;
375 ss << msgTag << " # max msg size : " << stringStats<global_size_t>(origComm, numActiveProcs, maxMsg) << std::endl;
376
377
378 outstr = ss.str();
379
380#ifdef HAVE_MPI
381 int strLength = outstr.size();
382 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
383 if (origComm->getRank() != root)
384 outstr.resize(strLength);
385 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
386#endif
387
388 return outstr;
389 }
390
391 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
392 std::string PerfUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CommPattern(const Matrix& A, const std::string& msgTag, RCP<const ParameterList> /* params */) {
393 if (!CheckMatrix(A))
394 return "";
395
396 std::ostringstream out;
397
398 RCP<const Teuchos::Comm<int> > comm = A.getRowMap()->getComm();
399 int myRank = comm->getRank();
400
401 out << msgTag << " " << myRank << ":";
402
403 RCP<const Import> importer = (A.getCrsGraph() != Teuchos::null ? A.getCrsGraph()->getImporter() : Teuchos::null);
404 if (importer.is_null()) {
405 out << std::endl;
406 return out.str();
407 }
408
409 ArrayView<const int> exportPIDs = importer->getExportPIDs();
410
411 if (exportPIDs.size()) {
412 // NodeTE: exportPIDs is sorted but not unique ( 1 1 1 2 2 3 4 4 4 )
413 int neigh = exportPIDs[0];
414 GO weight = 1;
415 for (int i = 1; i < exportPIDs.size(); i++) {
416 if (exportPIDs[i] != exportPIDs[i-1]) {
417 out << " " << neigh << "(" << weight << ")";
418
419 neigh = exportPIDs[i];
420 weight = 1;
421
422 } else {
423 weight += 1;
424 }
425 }
426 out << " " << neigh << "(" << weight << ")" << std::endl;
427 }
428
429 return out.str();
430 }
431
432 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
433 bool PerfUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CheckMatrix(const Matrix& A) {
434 // We can only print statistics for matrices that have a crs graph. A
435 // potential issue is regarding Xpetra::TpetraBlockCrsMatrix which has no
436 // CrsGraph. It is held as a private data member by Xpetra::CrsMatrix,
437 // which itself is an Xpetra::Matrix. So we check directly whether the
438 // request for the graph throws.
439 bool hasCrsGraph = true;
440 try {
441 A.getCrsGraph();
442
443 } catch (...) {
444 hasCrsGraph = false;
445 }
446
447 return hasCrsGraph;
448 }
449
450} //namespace MueLu
451
452#endif // MUELU_PERFUTILS_DEF_HPP
MueLu_PerfUtils_decl.hpp
MueLu_maxAll
#define MueLu_maxAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:65
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:61
MueLu_minAll
#define MueLu_minAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:63
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
MueLu::PerfUtils::CheckMatrix
static bool CheckMatrix(const Matrix &A)
Definition: MueLu_PerfUtils_def.hpp:433
MueLu::PerfUtils::PrintImporterInfo
static std::string PrintImporterInfo(RCP< const Import > importer, const std::string &msgTag)
Definition: MueLu_PerfUtils_def.hpp:329
MueLu::PerfUtils::PrintMatrixInfo
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Definition: MueLu_PerfUtils_def.hpp:127
MueLu::PerfUtils::CommPattern
static std::string CommPattern(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Definition: MueLu_PerfUtils_def.hpp:392
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::cmp_less
bool cmp_less(typename Map::value_type &v1, typename Map::value_type &v2)
Definition: MueLu_PerfUtils_def.hpp:122
MueLu::stringStats
std::string stringStats(const RCP< const Teuchos::Comm< int > > &comm, int numActiveProcs, const Type &v, RCP< ParameterList > paramList=Teuchos::null)
Definition: MueLu_PerfUtils_def.hpp:89
MueLu::calculateStats
void calculateStats(Type &minVal, Type &maxVal, double &avgVal, double &devVal, int &minProc, int &maxProc, const RCP< const Teuchos::Comm< int > > &comm, int numActiveProcs, const Type &v)
Definition: MueLu_PerfUtils_def.hpp:68