46 #ifndef MUELU_AGGREGATEQUALITYESTIMATEFACTORY_DEF_HPP
47 #define MUELU_AGGREGATEQUALITYESTIMATEFACTORY_DEF_HPP
60 #include "MueLu_Utilities.hpp"
66 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 Input(currentLevel,
"A");
75 Input(currentLevel,
"Aggregates");
76 Input(currentLevel,
"CoarseMap");
79 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
93 #undef SET_VALID_ENTRY
99 return validParamList;
102 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 RCP<Matrix> A = Get<RCP<Matrix>>(currentLevel,
"A");
107 RCP<Aggregates> aggregates = Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
109 RCP<const Map> map = Get<RCP<const Map>>(currentLevel,
"CoarseMap");
111 assert(!aggregates->AggregatesCrossProcessors());
113 std::string mode = pL.get<std::string>(
"aggregate qualities: mode");
114 GetOStream(
Statistics1) <<
"AggregateQuality: mode " << mode << std::endl;
117 if (mode ==
"eigenvalue" || mode ==
"both") {
119 ComputeAggregateQualities(A, aggregates, aggregate_qualities);
120 OutputAggQualities(currentLevel, aggregate_qualities);
122 if (mode ==
"size" || mode ==
"both") {
124 ComputeAggregateSizes(A, aggregates, aggregate_sizes);
125 Set(currentLevel,
"AggregateSizes", aggregate_sizes);
126 OutputAggSizes(currentLevel, aggregate_sizes);
128 Set(currentLevel,
"AggregateQualities", aggregate_qualities);
131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
147 aggsToIndices =
ArrayRCP<LO>(numAggs + LO_ONE, LO_ZERO);
149 for (
LO i = 0; i < numAggs; ++i) {
150 aggsToIndices[i + LO_ONE] = aggsToIndices[i] + aggSizes[i];
156 LO numNodes = vertex2AggId->getLocalLength();
158 std::vector<LO> vertexInsertionIndexByAgg(numNodes, LO_ZERO);
160 for (
LO i = 0; i < numNodes; ++i) {
161 LO aggId = vertex2AggIdData[i];
162 if (aggId < 0 || aggId >= numAggs)
continue;
164 aggSortedVertices[aggsToIndices[aggId] + vertexInsertionIndexByAgg[aggId]] = i;
165 vertexInsertionIndexByAgg[aggId]++;
169 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
172 const SC SCALAR_TWO = SCALAR_ONE + SCALAR_ONE;
185 std::string algostr = pL.
get<std::string>(
"aggregate qualities: algorithm");
186 MT zeroThreshold = Teuchos::as<MT>(pL.
get<
double>(
"aggregate qualities: zero threshold"));
187 enum AggAlgo { ALG_FORWARD = 0,
190 if (algostr ==
"forward") {
192 GetOStream(
Statistics1) <<
"AggregateQuality: Using 'forward' algorithm" << std::endl;
193 }
else if (algostr ==
"reverse") {
195 GetOStream(
Statistics1) <<
"AggregateQuality: Using 'reverse' algorithm" << std::endl;
200 bool check_symmetry = pL.
get<
bool>(
"aggregate qualities: check symmetry");
201 if (check_symmetry) {
203 x->Xpetra_randomize();
211 tmp->norm2(tmp_norm());
214 tmp->norm2(x_norm());
216 if (tmp_norm[0] > 1e-10 * x_norm[0]) {
217 std::string transpose_string =
"transpose";
221 assert(A->getMap()->isSameAs(*(AT->getMap())));
233 ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
234 ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
249 for (
LO aggId = LO_ZERO; aggId < numAggs; ++aggId) {
250 LO aggSize = aggSizes[aggId];
251 DenseMatrix A_aggPart(aggSize, aggSize,
true);
252 DenseVector offDiagonalAbsoluteSums(aggSize,
true);
255 for (
LO idx = LO_ZERO; idx < aggSize; ++idx) {
256 LO nodeId = aggSortedVertices[idx + aggsToIndices[aggId]];
257 A->getLocalRowView(nodeId, rowIndices, rowValues);
258 AT->getLocalRowView(nodeId, rowIndices, colValues);
261 for (
LO elem = LO_ZERO; elem < rowIndices.
size(); ++elem) {
262 LO nodeId2 = rowIndices[elem];
263 SC val = (rowValues[elem] + colValues[elem]) / SCALAR_TWO;
265 LO idxInAgg = -LO_ONE;
270 for (
LO idx2 = LO_ZERO; idx2 < aggSize; ++idx2) {
271 if (aggSortedVertices[idx2 + aggsToIndices[aggId]] == nodeId2) {
278 if (idxInAgg == -LO_ONE) {
291 DenseMatrix A_aggPartDiagonal(aggSize, aggSize,
true);
292 MT diag_sum = MT_ZERO;
293 for (
int i = 0; i < aggSize; ++i) {
298 DenseMatrix ones(aggSize, aggSize,
false);
299 ones.putScalar(MT_ONE);
303 DenseMatrix tmp(aggSize, aggSize,
false);
304 DenseMatrix topMatrix(A_aggPartDiagonal);
310 DenseMatrix bottomMatrix(A_aggPart);
311 MT matrixNorm = A_aggPart.normInf();
316 for (
int i = 0; i < aggSize; ++i) {
317 bottomMatrix(i, i) -= offDiagonalAbsoluteSums(i) + boost;
322 DenseVector alpha_real(aggSize,
false);
323 DenseVector alpha_imag(aggSize,
false);
324 DenseVector beta(aggSize,
false);
326 DenseVector workArray(14 * (aggSize + 1),
false);
335 const char compute_flag =
'N';
336 if (algo == ALG_FORWARD) {
338 myLapack.
GGES(compute_flag, compute_flag, compute_flag, ptr2func, aggSize,
339 topMatrix.values(), aggSize, bottomMatrix.values(), aggSize, &sdim,
340 alpha_real.values(), alpha_imag.values(), beta.values(), vl, aggSize,
341 vr, aggSize, workArray.values(), workArray.length(), bwork,
345 MT maxEigenVal = MT_ZERO;
346 for (
int i = LO_ZERO; i < aggSize; ++i) {
349 maxEigenVal = std::max(maxEigenVal, alpha_real[i] / beta[i]);
352 (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE + MT_ONE) * maxEigenVal;
356 myLapack.
GGES(compute_flag, compute_flag, compute_flag, ptr2func, aggSize,
357 bottomMatrix.values(), aggSize, topMatrix.values(), aggSize, &sdim,
358 alpha_real.values(), alpha_imag.values(), beta.values(), vl, aggSize,
359 vr, aggSize, workArray.values(), workArray.length(), bwork,
364 MT minEigenVal = MT_ZERO;
366 for (
int i = LO_ZERO; i < aggSize; ++i) {
367 MT ev = alpha_real[i] / beta[i];
368 if (ev > zeroThreshold) {
369 if (minEigenVal == MT_ZERO)
372 minEigenVal = std::min(minEigenVal, ev);
375 if (minEigenVal == MT_ZERO)
378 (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE + MT_ONE) / minEigenVal;
383 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
387 magnitudeType good_agg_thresh = Teuchos::as<magnitudeType>(pL.
get<
double>(
"aggregate qualities: good aggregate threshold"));
395 MT mean_bad_agg = 0.0;
396 MT mean_good_agg = 0.0;
398 for (
size_t i = 0; i < agg_qualities->getLocalLength(); ++i) {
399 if (data[i] > good_agg_thresh) {
401 mean_bad_agg += data[i];
403 mean_good_agg += data[i];
405 worst_agg = std::max(worst_agg, data[i]);
408 if (num_bad_aggs > 0) mean_bad_agg /= num_bad_aggs;
409 mean_good_agg /= agg_qualities->getLocalLength() - num_bad_aggs;
411 if (num_bad_aggs == 0) {
412 GetOStream(
Statistics1) <<
"All aggregates passed the quality measure. Worst aggregate had quality " << worst_agg <<
". Mean aggregate quality " << mean_good_agg <<
"." << std::endl;
414 GetOStream(
Statistics1) << num_bad_aggs <<
" out of " << agg_qualities->getLocalLength() <<
" did not pass the quality measure. Worst aggregate had quality " << worst_agg <<
". "
415 <<
"Mean bad aggregate quality " << mean_bad_agg <<
". Mean good aggregate quality " << mean_good_agg <<
"." << std::endl;
418 if (pL.
get<
bool>(
"aggregate qualities: file output")) {
419 std::string filename = pL.
get<std::string>(
"aggregate qualities: file base") +
"." + std::to_string(level.
GetLevelID());
424 const auto n = size_t(agg_qualities->getLocalLength());
429 for (
size_t i = 0; i < n; ++i) {
430 tmp.push_back(data[i]);
433 std::sort(tmp.begin(), tmp.end());
437 GetOStream(
Statistics1) <<
"AGG QUALITY HEADER : | LEVEL | TOTAL |";
438 for (
auto percent : percents) {
439 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) << 100.0 * percent <<
"% |";
444 for (
auto percent : percents) {
445 size_t i = size_t(n * percent);
446 i = i < n ? i : n - 1u;
448 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) << tmp[i] <<
" |";
454 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
456 ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
457 ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
460 auto data = agg_sizes->getDataNonConst(0);
461 for (
LO i = 0; i < (
LO)aggSizes.
size(); i++)
462 data[i] = aggSizes[i];
465 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
472 if (pL.
get<
bool>(
"aggregate qualities: file output")) {
473 std::string filename = pL.
get<std::string>(
"aggregate qualities: file base") +
".sizes." + std::to_string(level.
GetLevelID());
478 size_t n = (size_t)agg_sizes->getLocalLength();
483 for (
size_t i = 0; i < n; ++i) {
484 tmp.push_back(Teuchos::as<MT>(data[i]));
487 std::sort(tmp.begin(), tmp.end());
491 GetOStream(
Statistics1) <<
"AGG SIZE HEADER : | LEVEL | TOTAL |";
492 for (
auto percent : percents) {
493 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) << 100.0 * percent <<
"% |";
498 for (
auto percent : percents) {
499 size_t i = size_t(n * percent);
500 i = i < n ? i : n - 1u;
502 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) << tmp[i] <<
" |";
510 #endif // MUELU_AGGREGATEQUALITYESTIMATE_DEF_HPP
void OutputAggQualities(const Level &level, RCP< const Xpetra::MultiVector< magnitudeType, LO, GO, Node >> agg_qualities) const
virtual ~AggregateQualityEstimateFactory()
Destructor.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
KOKKOS_INLINE_FUNCTION LO GetNumAggregates() const
void ComputeAggregateSizes(RCP< const Matrix > A, RCP< const Aggregates > aggs, RCP< LocalOrdinalVector > agg_sizes) const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
static magnitudeType real(T a)
#define SET_VALID_ENTRY(name)
void DeclareInput(Level ¤tLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void GGES(const char &JOBVL, const char &JOBVR, const char &SORT, OrdinalType &(*ptr2func)(ScalarType *, ScalarType *, ScalarType *), const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, OrdinalType *sdim, MagnitudeType *ALPHAR, MagnitudeType *ALPHAI, MagnitudeType *BETA, ScalarType *VL, const OrdinalType &ldvl, ScalarType *VR, const OrdinalType &ldvr, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *BWORK, OrdinalType *info) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
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)
Teuchos::ArrayRCP< LocalOrdinal > ComputeAggregateSizesArrayRCP(bool forceRecompute=false) const
Compute sizes of aggregates.
void ComputeAggregateQualities(RCP< const Matrix > A, RCP< const Aggregates > aggs, RCP< Xpetra::MultiVector< magnitudeType, LO, GO, Node >> agg_qualities) const
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
void OutputAggSizes(const Level &level, RCP< const LocalOrdinalVector > agg_sizes) const
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
static magnitudeType magnitude(T a)
void Build(Level ¤tLevel) const
Build aggregate quality esimates with this factory.
static void ConvertAggregatesData(RCP< const Aggregates > aggs, ArrayRCP< LO > &aggSortedVertices, ArrayRCP< LO > &aggsToIndices, ArrayRCP< LO > &aggSizes)
Build aggregate quality esimates with this factory.
int GetLevelID() const
Return level number.
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
AggregateQualityEstimateFactory()
Constructor.