MueLu  Version of the Day
MueLu_AggregateQualityEstimateFactory_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_AGGREGATEQUALITYESTIMATEFACTORY_DEF_HPP
47 #define MUELU_AGGREGATEQUALITYESTIMATEFACTORY_DEF_HPP
48 #include <iomanip>
50 
51 #include "MueLu_Level.hpp"
52 
53 #include <Teuchos_SerialDenseVector.hpp>
54 #include <Teuchos_LAPACK.hpp>
55 
57 #include <Xpetra_IO.hpp>
58 #include "MueLu_MasterList.hpp"
59 #include "MueLu_Monitor.hpp"
60 #include "MueLu_FactoryManager.hpp"
61 #include "MueLu_Utilities.hpp"
62 
63 #include <vector>
64 
65 namespace MueLu {
66 
67  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  { }
70 
71  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
73 
74  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
76 
77  Input(currentLevel, "A");
78  Input(currentLevel, "Aggregates");
79  Input(currentLevel, "CoarseMap");
80 
81  }
82 
83  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
85  RCP<ParameterList> validParamList = rcp(new ParameterList());
86 
87 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
88  SET_VALID_ENTRY("aggregate qualities: good aggregate threshold");
89  SET_VALID_ENTRY("aggregate qualities: file output");
90  SET_VALID_ENTRY("aggregate qualities: file base");
91  SET_VALID_ENTRY("aggregate qualities: check symmetry");
92  SET_VALID_ENTRY("aggregate qualities: algorithm");
93  SET_VALID_ENTRY("aggregate qualities: zero threshold");
94  SET_VALID_ENTRY("aggregate qualities: percentiles");
95  SET_VALID_ENTRY("aggregate qualities: mode");
96 
97 #undef SET_VALID_ENTRY
98 
99  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
100  validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Generating factory of the aggregates");
101  validParamList->set< RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
102 
103  return validParamList;
104  }
105 
106 
107  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
109 
110  FactoryMonitor m(*this, "Build", currentLevel);
111 
112  RCP<Matrix> A = Get<RCP<Matrix>>(currentLevel, "A");
113  RCP<Aggregates> aggregates = Get<RCP<Aggregates>>(currentLevel, "Aggregates");
114 
115  RCP<const Map> map = Get< RCP<const Map> >(currentLevel, "CoarseMap");
116 
117 
118  assert(!aggregates->AggregatesCrossProcessors());
119  ParameterList pL = GetParameterList();
120  std::string mode = pL.get<std::string>("aggregate qualities: mode");
121  GetOStream(Statistics1) << "AggregateQuality: mode "<<mode << std::endl;
122 
123  RCP<Xpetra::MultiVector<magnitudeType,LO,GO,NO>> aggregate_qualities;
124  if(mode == "eigenvalue" || mode == "both") {
125  aggregate_qualities = Xpetra::MultiVectorFactory<magnitudeType,LO,GO,NO>::Build(map, 1);
126  ComputeAggregateQualities(A, aggregates, aggregate_qualities);
127  OutputAggQualities(currentLevel, aggregate_qualities);
128  }
129  if(mode == "size" || mode =="both") {
130  RCP<LocalOrdinalVector> aggregate_sizes = Xpetra::VectorFactory<LO,LO,GO,NO>::Build(map);
131  ComputeAggregateSizes(A,aggregates,aggregate_sizes);
132  Set(currentLevel, "AggregateSizes",aggregate_sizes);
133  OutputAggSizes(currentLevel, aggregate_sizes);
134  }
135  Set(currentLevel, "AggregateQualities", aggregate_qualities);
136 
137 
138  }
139 
140  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141  void AggregateQualityEstimateFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ConvertAggregatesData(RCP<const Aggregates> aggs, ArrayRCP<LO>& aggSortedVertices, ArrayRCP<LO>& aggsToIndices, ArrayRCP<LO>& aggSizes) {
142 
143  // Reorder local aggregate information into a format amenable to computing
144  // per-aggregate quantities. Specifically, we compute a format
145  // similar to compressed sparse row format for sparse matrices in which
146  // we store all the local vertices in a single array in blocks corresponding
147  // to aggregates. (This array is aggSortedVertices.) We then store a second
148  // array (aggsToIndices) whose k-th element stores the index of the first
149  // vertex in aggregate k in the array aggSortedVertices.
150 
151  const LO LO_ZERO = Teuchos::OrdinalTraits<LO>::zero();
152  const LO LO_ONE = Teuchos::OrdinalTraits<LO>::one();
153 
154  LO numAggs = aggs->GetNumAggregates();
155  aggSizes = aggs->ComputeAggregateSizes();
156 
157  aggsToIndices = ArrayRCP<LO>(numAggs+LO_ONE,LO_ZERO);
158 
159  for (LO i=0;i<numAggs;++i) {
160  aggsToIndices[i+LO_ONE] = aggsToIndices[i] + aggSizes[i];
161  }
162 
163  const RCP<LOMultiVector> vertex2AggId = aggs->GetVertex2AggId();
164  const ArrayRCP<const LO> vertex2AggIdData = vertex2AggId->getData(0);
165 
166  LO numNodes = vertex2AggId->getLocalLength();
167  aggSortedVertices = ArrayRCP<LO>(numNodes,-LO_ONE);
168  std::vector<LO> vertexInsertionIndexByAgg(numNodes,LO_ZERO);
169 
170  for (LO i=0;i<numNodes;++i) {
171 
172  LO aggId = vertex2AggIdData[i];
173  if (aggId<0 || aggId>=numAggs) continue;
174 
175  aggSortedVertices[aggsToIndices[aggId]+vertexInsertionIndexByAgg[aggId]] = i;
176  vertexInsertionIndexByAgg[aggId]++;
177 
178  }
179 
180 
181  }
182 
183  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
184  void AggregateQualityEstimateFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeAggregateQualities(RCP<const Matrix> A, RCP<const Aggregates> aggs, RCP<Xpetra::MultiVector<magnitudeType,LO,GO,Node>> agg_qualities) const {
185 
186  const SC SCALAR_ONE = Teuchos::ScalarTraits<SC>::one();
187  const SC SCALAR_TWO = SCALAR_ONE + SCALAR_ONE;
188 
189  const LO LO_ZERO = Teuchos::OrdinalTraits<LO>::zero();
190  const LO LO_ONE = Teuchos::OrdinalTraits<LO>::one();
191 
192  using MT = magnitudeType;
193  const MT MT_ZERO = Teuchos::ScalarTraits<MT>::zero();
194  const MT MT_ONE = Teuchos::ScalarTraits<MT>::one();
195  ParameterList pL = GetParameterList();
196 
197  RCP<const Matrix> AT = A;
198 
199  // Algorithm check
200  std::string algostr = pL.get<std::string>("aggregate qualities: algorithm");
201  MT zeroThreshold = Teuchos::as<MT>(pL.get<double>("aggregate qualities: zero threshold"));
202  enum AggAlgo {ALG_FORWARD=0, ALG_REVERSE};
203  AggAlgo algo;
204  if(algostr == "forward") {algo = ALG_FORWARD; GetOStream(Statistics1) << "AggregateQuality: Using 'forward' algorithm" << std::endl;}
205  else if(algostr == "reverse") {algo = ALG_REVERSE; GetOStream(Statistics1) << "AggregateQuality: Using 'reverse' algorithm" << std::endl;}
206  else {
207  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "\"algorithm\" must be one of (forward|reverse)");
208  }
209 
210  bool check_symmetry = pL.get<bool>("aggregate qualities: check symmetry");
211  if (check_symmetry) {
212 
213  RCP<MultiVector> x = MultiVectorFactory::Build(A->getMap(), 1, false);
214  x->Xpetra_randomize();
215 
216  RCP<MultiVector> tmp = MultiVectorFactory::Build(A->getMap(), 1, false);
217 
218  A->apply(*x, *tmp, Teuchos::NO_TRANS); // tmp now stores A*x
219  A->apply(*x, *tmp, Teuchos::TRANS, -SCALAR_ONE, SCALAR_ONE); // tmp now stores A*x - A^T*x
220 
221  Array<magnitudeType> tmp_norm(1);
222  tmp->norm2(tmp_norm());
223 
224  Array<magnitudeType> x_norm(1);
225  tmp->norm2(x_norm());
226 
227  if (tmp_norm[0] > 1e-10*x_norm[0]) {
228  std::string transpose_string = "transpose";
229  RCP<ParameterList> whatever;
230  AT = Utilities::Transpose(*rcp_const_cast<Matrix>(A), true, transpose_string, whatever);
231 
232  assert(A->getMap()->isSameAs( *(AT->getMap()) ));
233  }
234 
235  }
236 
237  // Reorder local aggregate information into a format amenable to computing
238  // per-aggregate quantities. Specifically, we compute a format
239  // similar to compressed sparse row format for sparse matrices in which
240  // we store all the local vertices in a single array in blocks corresponding
241  // to aggregates. (This array is aggSortedVertices.) We then store a second
242  // array (aggsToIndices) whose k-th element stores the index of the first
243  // vertex in aggregate k in the array aggSortedVertices.
244 
245  ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
246  ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
247 
248  LO numAggs = aggs->GetNumAggregates();
249 
250  // Compute the per-aggregate quality estimate
251 
252  typedef Teuchos::SerialDenseMatrix<LO,MT> DenseMatrix;
253  typedef Teuchos::SerialDenseVector<LO,MT> DenseVector;
254 
255  ArrayView<const LO> rowIndices;
256  ArrayView<const SC> rowValues;
257  ArrayView<const SC> colValues;
258  Teuchos::LAPACK<LO,MT> myLapack;
259 
260  // Iterate over each aggregate to compute the quality estimate
261  for (LO aggId=LO_ZERO; aggId<numAggs; ++aggId) {
262 
263  LO aggSize = aggSizes[aggId];
264  DenseMatrix A_aggPart(aggSize, aggSize, true);
265  DenseVector offDiagonalAbsoluteSums(aggSize, true);
266 
267  // Iterate over each node in the aggregate
268  for (LO idx=LO_ZERO; idx<aggSize; ++idx) {
269 
270  LO nodeId = aggSortedVertices[idx+aggsToIndices[aggId]];
271  A->getLocalRowView(nodeId, rowIndices, rowValues);
272  AT->getLocalRowView(nodeId, rowIndices, colValues);
273 
274  // Iterate over each element in the row corresponding to the current node
275  for (LO elem=LO_ZERO; elem<rowIndices.size();++elem) {
276 
277  LO nodeId2 = rowIndices[elem];
278  SC val = (rowValues[elem]+colValues[elem])/SCALAR_TWO;
279 
280  LO idxInAgg = -LO_ONE; // -1 if element is not in aggregate
281 
282  // Check whether the element belongs in the aggregate. If it does
283  // find, its index. Otherwise, add it's value to the off diagonal
284  // sums
285  for (LO idx2=LO_ZERO; idx2<aggSize; ++idx2) {
286 
287  if (aggSortedVertices[idx2+aggsToIndices[aggId]] == nodeId2) {
288 
289  // Element does belong to aggregate
290  idxInAgg = idx2;
291  break;
292 
293  }
294 
295  }
296 
297  if (idxInAgg == -LO_ONE) { // Element does not belong to aggregate
298 
299  offDiagonalAbsoluteSums[idx] += Teuchos::ScalarTraits<SC>::magnitude(val);
300 
301  } else { // Element does belong to aggregate
302 
303  A_aggPart(idx,idxInAgg) = Teuchos::ScalarTraits<SC>::real(val);
304 
305  }
306 
307  }
308 
309  }
310 
311  // Construct a diagonal matrix consisting of the diagonal
312  // of A_aggPart
313  DenseMatrix A_aggPartDiagonal(aggSize, aggSize, true);
314  MT diag_sum = MT_ZERO;
315  for (int i=0;i<aggSize;++i) {
316  A_aggPartDiagonal(i,i) = Teuchos::ScalarTraits<SC>::real(A_aggPart(i,i));
317  diag_sum += Teuchos::ScalarTraits<SC>::real(A_aggPart(i,i));
318  }
319 
320  DenseMatrix ones(aggSize, aggSize, false);
321  ones.putScalar(MT_ONE);
322 
323  // Compute matrix on top of generalized Rayleigh quotient
324  // topMatrix = A_aggPartDiagonal - A_aggPartDiagonal*ones*A_aggPartDiagonal/diag_sum;
325  DenseMatrix tmp(aggSize, aggSize, false);
326  DenseMatrix topMatrix(A_aggPartDiagonal);
327 
328  tmp.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, MT_ONE, ones, A_aggPartDiagonal, MT_ZERO);
329  topMatrix.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -MT_ONE/diag_sum, A_aggPartDiagonal, tmp, MT_ONE);
330 
331  // Compute matrix on bottom of generalized Rayleigh quotient
332  DenseMatrix bottomMatrix(A_aggPart);
333  MT matrixNorm = A_aggPart.normInf();
334 
335  // Forward mode: Include a small perturbation to the bottom matrix to make it nonsingular
336  const MT boost = (algo == ALG_FORWARD) ? (-1e4*Teuchos::ScalarTraits<MT>::eps()*matrixNorm) : MT_ZERO;
337 
338  for (int i=0;i<aggSize;++i){
339  bottomMatrix(i,i) -= offDiagonalAbsoluteSums(i) + boost;
340  }
341 
342  // Compute generalized eigenvalues
343  LO sdim, info;
344  DenseVector alpha_real(aggSize, false);
345  DenseVector alpha_imag(aggSize, false);
346  DenseVector beta(aggSize, false);
347 
348  DenseVector workArray(14*(aggSize+1), false);
349 
350  LO (*ptr2func)(MT*, MT*, MT*);
351  ptr2func = NULL;
352  LO* bwork = NULL;
353  MT* vl = NULL;
354  MT* vr = NULL;
355 
356  const char compute_flag ='N';
357  if(algo == ALG_FORWARD) {
358  // Forward: Solve the generalized eigenvalue problem as is
359  myLapack.GGES(compute_flag,compute_flag,compute_flag,ptr2func,aggSize,
360  topMatrix.values(),aggSize,bottomMatrix.values(),aggSize,&sdim,
361  alpha_real.values(),alpha_imag.values(),beta.values(),vl,aggSize,
362  vr,aggSize,workArray.values(),workArray.length(),bwork,
363  &info);
364  TEUCHOS_ASSERT(info == LO_ZERO);
365 
366  MT maxEigenVal = MT_ZERO;
367  for (int i=LO_ZERO;i<aggSize;++i) {
368  // NOTE: In theory, the eigenvalues should be nearly real
369  //TEUCHOS_ASSERT(fabs(alpha_imag[i]) <= 1e-8*fabs(alpha_real[i])); // Eigenvalues should be nearly real
370  maxEigenVal = std::max(maxEigenVal, alpha_real[i]/beta[i]);
371  }
372 
373  (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE+MT_ONE)*maxEigenVal;
374  }
375  else {
376  // Reverse: Swap the top and bottom matrices for the generalized eigenvalue problem
377  // This is trickier, since we need to grab the smallest non-zero eigenvalue and invert it.
378  myLapack.GGES(compute_flag,compute_flag,compute_flag,ptr2func,aggSize,
379  bottomMatrix.values(),aggSize,topMatrix.values(),aggSize,&sdim,
380  alpha_real.values(),alpha_imag.values(),beta.values(),vl,aggSize,
381  vr,aggSize,workArray.values(),workArray.length(),bwork,
382  &info);
383 
384  TEUCHOS_ASSERT(info == LO_ZERO);
385 
386  MT minEigenVal = MT_ZERO;
387 
388  for (int i=LO_ZERO;i<aggSize;++i) {
389  MT ev = alpha_real[i] / beta[i];
390  if(ev > zeroThreshold) {
391  if (minEigenVal == MT_ZERO) minEigenVal = ev;
392  else minEigenVal = std::min(minEigenVal,ev);
393  }
394  }
395  if(minEigenVal == MT_ZERO) (agg_qualities->getDataNonConst(0))[aggId] = Teuchos::ScalarTraits<MT>::rmax();
396  else (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE+MT_ONE) / minEigenVal;
397  }
398  }//end aggId loop
399  }
400 
401  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
402  void AggregateQualityEstimateFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::OutputAggQualities(const Level& level, RCP<const Xpetra::MultiVector<magnitudeType,LO,GO,Node>> agg_qualities) const {
403 
404  ParameterList pL = GetParameterList();
405 
406  magnitudeType good_agg_thresh = Teuchos::as<magnitudeType>(pL.get<double>("aggregate qualities: good aggregate threshold"));
407  using MT = magnitudeType;
408 
409  ArrayRCP<const MT> data = agg_qualities->getData(0);
410 
411  LO num_bad_aggs = 0;
412  MT worst_agg = 0.0;
413 
414  MT mean_bad_agg = 0.0;
415  MT mean_good_agg = 0.0;
416 
417 
418  for (size_t i=0;i<agg_qualities->getLocalLength();++i) {
419 
420  if (data[i] > good_agg_thresh) {
421  num_bad_aggs++;
422  mean_bad_agg += data[i];
423  }
424  else {
425  mean_good_agg += data[i];
426  }
427  worst_agg = std::max(worst_agg, data[i]);
428  }
429 
430 
431  if (num_bad_aggs > 0) mean_bad_agg /= num_bad_aggs;
432  mean_good_agg /= agg_qualities->getLocalLength() - num_bad_aggs;
433 
434  if (num_bad_aggs == 0) {
435  GetOStream(Statistics1) << "All aggregates passed the quality measure. Worst aggregate had quality " << worst_agg << ". Mean aggregate quality " << mean_good_agg << "." << std::endl;
436  } else {
437  GetOStream(Statistics1) << num_bad_aggs << " out of " << agg_qualities->getLocalLength() << " did not pass the quality measure. Worst aggregate had quality " << worst_agg << ". "
438  << "Mean bad aggregate quality " << mean_bad_agg << ". Mean good aggregate quality " << mean_good_agg << "." << std::endl;
439  }
440 
441  if (pL.get<bool>("aggregate qualities: file output")) {
442  std::string filename = pL.get<std::string>("aggregate qualities: file base")+"."+std::to_string(level.GetLevelID());
443  Xpetra::IO<magnitudeType,LO,GO,Node>::Write(filename, *agg_qualities);
444  }
445 
446  {
447  const auto n = size_t(agg_qualities->getLocalLength());
448 
449  std::vector<MT> tmp;
450  tmp.reserve(n);
451 
452  for (size_t i=0; i<n; ++i) {
453  tmp.push_back(data[i]);
454  }
455 
456  std::sort(tmp.begin(), tmp.end());
457 
458  Teuchos::ArrayView<const double> percents = pL.get<Teuchos::Array<double> >("aggregate qualities: percentiles")();
459 
460  GetOStream(Statistics1) << "AGG QUALITY HEADER : | LEVEL | TOTAL |";
461  for (auto percent : percents) {
462  GetOStream(Statistics1) << std::fixed << std::setprecision(4) <<100.0*percent << "% |";
463  }
464  GetOStream(Statistics1) << std::endl;
465 
466  GetOStream(Statistics1) << "AGG QUALITY PERCENTILES: | " << level.GetLevelID() << " | " << n << "|";
467  for (auto percent : percents) {
468  size_t i = size_t(n*percent);
469  i = i < n ? i : n-1u;
470  i = i > 0u ? i : 0u;
471  GetOStream(Statistics1) << std::fixed <<std::setprecision(4) << tmp[i] << " |";
472  }
473  GetOStream(Statistics1) << std::endl;
474 
475  }
476  }
477 
478 
479 
480 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
481  void AggregateQualityEstimateFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeAggregateSizes(RCP<const Matrix> A, RCP<const Aggregates> aggs, RCP<LocalOrdinalVector> agg_sizes) const {
482 
483  ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
484  ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
485 
486  // Iterate over each node in the aggregate
487  auto data = agg_sizes->getDataNonConst(0);
488  for (LO i=0; i<(LO)aggSizes.size(); i++)
489  data[i] = aggSizes[i];
490 }
491 
492 
493 
494 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
495  void AggregateQualityEstimateFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::OutputAggSizes(const Level& level, RCP<const LocalOrdinalVector> agg_sizes) const {
496 
497  ParameterList pL = GetParameterList();
498  using MT = magnitudeType;
499 
500  ArrayRCP<const LO> data = agg_sizes->getData(0);
501 
502 
503  if (pL.get<bool>("aggregate qualities: file output")) {
504  std::string filename = pL.get<std::string>("aggregate qualities: file base")+".sizes."+std::to_string(level.GetLevelID());
505  Xpetra::IO<LO,LO,GO,Node>::Write(filename, *agg_sizes);
506  }
507 
508  {
509  size_t n = (size_t)agg_sizes->getLocalLength();
510 
511  std::vector<MT> tmp;
512  tmp.reserve(n);
513 
514  for (size_t i=0; i<n; ++i) {
515  tmp.push_back(Teuchos::as<MT>(data[i]));
516  }
517 
518  std::sort(tmp.begin(), tmp.end());
519 
520  Teuchos::ArrayView<const double> percents = pL.get<Teuchos::Array<double> >("aggregate qualities: percentiles")();
521 
522  GetOStream(Statistics1) << "AGG SIZE HEADER : | LEVEL | TOTAL |";
523  for (auto percent : percents) {
524  GetOStream(Statistics1) << std::fixed << std::setprecision(4) <<100.0*percent << "% |";
525  }
526  GetOStream(Statistics1) << std::endl;
527 
528  GetOStream(Statistics1) << "AGG SIZE PERCENTILES: | " << level.GetLevelID() << " | " << n << "|";
529  for (auto percent : percents) {
530  size_t i = size_t(n*percent);
531  i = i < n ? i : n-1u;
532  i = i > 0u ? i : 0u;
533  GetOStream(Statistics1) << std::fixed <<std::setprecision(4) << tmp[i] << " |";
534  }
535  GetOStream(Statistics1) << std::endl;
536 
537  }
538  }
539 
540 
541 
542 } // namespace MueLu
543 
544 #endif // MUELU_AGGREGATEQUALITYESTIMATE_DEF_HPP
void OutputAggQualities(const Level &level, RCP< const Xpetra::MultiVector< magnitudeType, LO, GO, Node >> agg_qualities) const
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
Namespace for MueLu classes and methods.
#define SET_VALID_ENTRY(name)
void OutputAggSizes(const Level &level, RCP< const LocalOrdinalVector > agg_sizes) const
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void ComputeAggregateSizes(RCP< const Matrix > A, RCP< const Aggregates > aggs, RCP< LocalOrdinalVector > agg_sizes) const
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
void ComputeAggregateQualities(RCP< const Matrix > A, RCP< const Aggregates > aggs, RCP< Xpetra::MultiVector< magnitudeType, LO, GO, Node >> agg_qualities) const
void Build(Level &currentLevel) const
Build aggregate quality esimates with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static void ConvertAggregatesData(RCP< const Aggregates > aggs, ArrayRCP< LO > &aggSortedVertices, ArrayRCP< LO > &aggsToIndices, ArrayRCP< LO > &aggSizes)
Build aggregate quality esimates with this factory.
Exception throws to report errors in the internal logical of the program.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType