46#ifndef MUELU_AGGREGATEQUALITYESTIMATEFACTORY_DEF_HPP
47#define MUELU_AGGREGATEQUALITYESTIMATEFACTORY_DEF_HPP
53#include <Teuchos_SerialDenseVector.hpp>
54#include <Teuchos_LAPACK.hpp>
57#include <Xpetra_IO.hpp>
60#include "MueLu_FactoryManager.hpp"
61#include "MueLu_Utilities.hpp"
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 Input(currentLevel,
"A");
78 Input(currentLevel,
"Aggregates");
79 Input(currentLevel,
"CoarseMap");
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 RCP<ParameterList> validParamList = rcp(
new ParameterList());
87#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
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");
103 return validParamList;
107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 RCP<Matrix> A = Get<RCP<Matrix>>(currentLevel,
"A");
113 RCP<Aggregates> aggregates = Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
115 RCP<const Map> map = Get< RCP<const Map> >(currentLevel,
"CoarseMap");
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;
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);
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);
135 Set(currentLevel,
"AggregateQualities", aggregate_qualities);
140 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
151 const LO LO_ZERO = Teuchos::OrdinalTraits<LO>::zero();
152 const LO LO_ONE = Teuchos::OrdinalTraits<LO>::one();
154 LO numAggs = aggs->GetNumAggregates();
155 aggSizes = aggs->ComputeAggregateSizes();
157 aggsToIndices = ArrayRCP<LO>(numAggs+LO_ONE,LO_ZERO);
159 for (LO i=0;i<numAggs;++i) {
160 aggsToIndices[i+LO_ONE] = aggsToIndices[i] + aggSizes[i];
163 const RCP<LOMultiVector> vertex2AggId = aggs->GetVertex2AggId();
164 const ArrayRCP<const LO> vertex2AggIdData = vertex2AggId->getData(0);
166 LO numNodes = vertex2AggId->getLocalLength();
167 aggSortedVertices = ArrayRCP<LO>(numNodes,-LO_ONE);
168 std::vector<LO> vertexInsertionIndexByAgg(numNodes,LO_ZERO);
170 for (LO i=0;i<numNodes;++i) {
172 LO aggId = vertex2AggIdData[i];
173 if (aggId<0 || aggId>=numAggs)
continue;
175 aggSortedVertices[aggsToIndices[aggId]+vertexInsertionIndexByAgg[aggId]] = i;
176 vertexInsertionIndexByAgg[aggId]++;
183 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
186 const SC SCALAR_ONE = Teuchos::ScalarTraits<SC>::one();
187 const SC SCALAR_TWO = SCALAR_ONE + SCALAR_ONE;
189 const LO LO_ZERO = Teuchos::OrdinalTraits<LO>::zero();
190 const LO LO_ONE = Teuchos::OrdinalTraits<LO>::one();
193 const MT MT_ZERO = Teuchos::ScalarTraits<MT>::zero();
194 const MT MT_ONE = Teuchos::ScalarTraits<MT>::one();
195 ParameterList pL = GetParameterList();
197 RCP<const Matrix> AT = A;
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};
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;}
210 bool check_symmetry = pL.get<
bool>(
"aggregate qualities: check symmetry");
211 if (check_symmetry) {
213 RCP<MultiVector> x = MultiVectorFactory::Build(A->getMap(), 1,
false);
214 x->Xpetra_randomize();
216 RCP<MultiVector> tmp = MultiVectorFactory::Build(A->getMap(), 1,
false);
218 A->apply(*x, *tmp, Teuchos::NO_TRANS);
219 A->apply(*x, *tmp, Teuchos::TRANS, -SCALAR_ONE, SCALAR_ONE);
221 Array<magnitudeType> tmp_norm(1);
222 tmp->norm2(tmp_norm());
224 Array<magnitudeType> x_norm(1);
225 tmp->norm2(x_norm());
227 if (tmp_norm[0] > 1e-10*x_norm[0]) {
228 std::string transpose_string =
"transpose";
229 RCP<ParameterList> whatever;
232 assert(A->getMap()->isSameAs( *(AT->getMap()) ));
245 ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
246 ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
248 LO numAggs = aggs->GetNumAggregates();
252 typedef Teuchos::SerialDenseMatrix<LO,MT> DenseMatrix;
253 typedef Teuchos::SerialDenseVector<LO,MT> DenseVector;
255 ArrayView<const LO> rowIndices;
256 ArrayView<const SC> rowValues;
257 ArrayView<const SC> colValues;
258 Teuchos::LAPACK<LO,MT> myLapack;
261 for (LO aggId=LO_ZERO; aggId<numAggs; ++aggId) {
263 LO aggSize = aggSizes[aggId];
264 DenseMatrix A_aggPart(aggSize, aggSize,
true);
265 DenseVector offDiagonalAbsoluteSums(aggSize,
true);
268 for (LO idx=LO_ZERO; idx<aggSize; ++idx) {
270 LO nodeId = aggSortedVertices[idx+aggsToIndices[aggId]];
271 A->getLocalRowView(nodeId, rowIndices, rowValues);
272 AT->getLocalRowView(nodeId, rowIndices, colValues);
275 for (LO elem=LO_ZERO; elem<rowIndices.size();++elem) {
277 LO nodeId2 = rowIndices[elem];
278 SC val = (rowValues[elem]+colValues[elem])/SCALAR_TWO;
280 LO idxInAgg = -LO_ONE;
285 for (LO idx2=LO_ZERO; idx2<aggSize; ++idx2) {
287 if (aggSortedVertices[idx2+aggsToIndices[aggId]] == nodeId2) {
297 if (idxInAgg == -LO_ONE) {
299 offDiagonalAbsoluteSums[idx] += Teuchos::ScalarTraits<SC>::magnitude(val);
303 A_aggPart(idx,idxInAgg) = Teuchos::ScalarTraits<SC>::real(val);
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));
320 DenseMatrix ones(aggSize, aggSize,
false);
321 ones.putScalar(MT_ONE);
325 DenseMatrix tmp(aggSize, aggSize,
false);
326 DenseMatrix topMatrix(A_aggPartDiagonal);
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);
332 DenseMatrix bottomMatrix(A_aggPart);
333 MT matrixNorm = A_aggPart.normInf();
336 const MT boost = (algo == ALG_FORWARD) ? (-1e4*Teuchos::ScalarTraits<MT>::eps()*matrixNorm) : MT_ZERO;
338 for (
int i=0;i<aggSize;++i){
339 bottomMatrix(i,i) -= offDiagonalAbsoluteSums(i) + boost;
344 DenseVector alpha_real(aggSize,
false);
345 DenseVector alpha_imag(aggSize,
false);
346 DenseVector beta(aggSize,
false);
348 DenseVector workArray(14*(aggSize+1),
false);
350 LO (*ptr2func)(MT*, MT*, MT*);
356 const char compute_flag =
'N';
357 if(algo == ALG_FORWARD) {
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,
364 TEUCHOS_ASSERT(info == LO_ZERO);
366 MT maxEigenVal = MT_ZERO;
367 for (
int i=LO_ZERO;i<aggSize;++i) {
370 maxEigenVal = std::max(maxEigenVal, alpha_real[i]/beta[i]);
373 (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE+MT_ONE)*maxEigenVal;
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,
384 TEUCHOS_ASSERT(info == LO_ZERO);
386 MT minEigenVal = MT_ZERO;
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);
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;
401 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
404 ParameterList pL = GetParameterList();
406 magnitudeType good_agg_thresh = Teuchos::as<magnitudeType>(pL.get<
double>(
"aggregate qualities: good aggregate threshold"));
409 ArrayRCP<const MT> data = agg_qualities->getData(0);
414 MT mean_bad_agg = 0.0;
415 MT mean_good_agg = 0.0;
418 for (
size_t i=0;i<agg_qualities->getLocalLength();++i) {
420 if (data[i] > good_agg_thresh) {
422 mean_bad_agg += data[i];
425 mean_good_agg += data[i];
427 worst_agg = std::max(worst_agg, data[i]);
431 if (num_bad_aggs > 0) mean_bad_agg /= num_bad_aggs;
432 mean_good_agg /= agg_qualities->getLocalLength() - num_bad_aggs;
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;
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;
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);
447 const auto n = size_t(agg_qualities->getLocalLength());
452 for (
size_t i=0; i<n; ++i) {
453 tmp.push_back(data[i]);
456 std::sort(tmp.begin(), tmp.end());
458 Teuchos::ArrayView<const double> percents = pL.get<Teuchos::Array<double> >(
"aggregate qualities: percentiles")();
460 GetOStream(
Statistics1) <<
"AGG QUALITY HEADER : | LEVEL | TOTAL |";
461 for (
auto percent : percents) {
462 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) <<100.0*percent <<
"% |";
467 for (
auto percent : percents) {
468 size_t i = size_t(n*percent);
469 i = i < n ? i : n-1u;
471 GetOStream(
Statistics1) << std::fixed <<std::setprecision(4) << tmp[i] <<
" |";
480template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
483 ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
484 ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
487 auto data = agg_sizes->getDataNonConst(0);
488 for (LO i=0; i<(LO)aggSizes.size(); i++)
489 data[i] = aggSizes[i];
494template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
497 ParameterList pL = GetParameterList();
500 ArrayRCP<const LO> data = agg_sizes->getData(0);
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);
509 size_t n = (size_t)agg_sizes->getLocalLength();
514 for (
size_t i=0; i<n; ++i) {
515 tmp.push_back(Teuchos::as<MT>(data[i]));
518 std::sort(tmp.begin(), tmp.end());
520 Teuchos::ArrayView<const double> percents = pL.get<Teuchos::Array<double> >(
"aggregate qualities: percentiles")();
522 GetOStream(
Statistics1) <<
"AGG SIZE HEADER : | LEVEL | TOTAL |";
523 for (
auto percent : percents) {
524 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) <<100.0*percent <<
"% |";
529 for (
auto percent : percents) {
530 size_t i = size_t(n*percent);
531 i = i < n ? i : n-1u;
533 GetOStream(
Statistics1) << std::fixed <<std::setprecision(4) << tmp[i] <<
" |";
#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 ¤tLevel) 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 ¤tLevel) 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
virtual ~AggregateQualityEstimateFactory()
Destructor.
AggregateQualityEstimateFactory()
Constructor.
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 > ¶ms=Teuchos::null)
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.