MueLu Version of the Day
Loading...
Searching...
No Matches
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
65namespace MueLu {
66
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
480template <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
494template <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
#define SET_VALID_ENTRY(name)
void OutputAggQualities(const Level &level, RCP< const Xpetra::MultiVector< magnitudeType, LO, GO, Node > > agg_qualities) const
void OutputAggSizes(const Level &level, RCP< const LocalOrdinalVector > agg_sizes) 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.
void ComputeAggregateQualities(RCP< const Matrix > A, RCP< const Aggregates > aggs, RCP< Xpetra::MultiVector< magnitudeType, LO, GO, Node > > agg_qualities) const
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
static void ConvertAggregatesData(RCP< const Aggregates > aggs, ArrayRCP< LO > &aggSortedVertices, ArrayRCP< LO > &aggsToIndices, ArrayRCP< LO > &aggSizes)
Build aggregate quality esimates with this factory.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
void ComputeAggregateSizes(RCP< const Matrix > A, RCP< const Aggregates > aggs, RCP< LocalOrdinalVector > agg_sizes) const
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
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)
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.