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)) 97 #undef SET_VALID_ENTRY 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] <<
" |";
480 template <
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];
494 template <
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] <<
" |";
544 #endif // MUELU_AGGREGATEQUALITYESTIMATE_DEF_HPP void OutputAggQualities(const Level &level, RCP< const Xpetra::MultiVector< magnitudeType, LO, GO, Node >> agg_qualities) const
virtual ~AggregateQualityEstimateFactory()
Destructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
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 ¤tLevel) const
Specifies the data that this class needs, and the factories that generate that data.
int GetLevelID() const
Return level number.
Class that holds all level-specific information.
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 > ¶ms=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 ¤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.
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
AggregateQualityEstimateFactory()
Constructor.