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)) 95 #undef SET_VALID_ENTRY 97 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
98 validParamList->set< RCP<const FactoryBase> >(
"Aggregates", Teuchos::null,
"Generating factory of the aggregates");
99 validParamList->set< RCP<const FactoryBase> >(
"CoarseMap", Teuchos::null,
"Generating factory of the coarse map");
101 return validParamList;
105 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 RCP<Matrix> A = Get<RCP<Matrix>>(currentLevel,
"A");
111 RCP<Aggregates> aggregates = Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
113 RCP<const Map> map = Get< RCP<const Map> >(currentLevel,
"CoarseMap");
114 RCP<Xpetra::MultiVector<magnitudeType,LO,GO,NO>> aggregate_qualities = Xpetra::MultiVectorFactory<magnitudeType,LO,GO,NO>::Build(map, 1);
116 assert(!aggregates->AggregatesCrossProcessors());
117 ComputeAggregateQualities(A, aggregates, aggregate_qualities);
119 Set(currentLevel,
"AggregateQualities", aggregate_qualities);
121 OutputAggQualities(currentLevel, aggregate_qualities);
125 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
136 const LO LO_ZERO = Teuchos::OrdinalTraits<LO>::zero();
137 const LO LO_ONE = Teuchos::OrdinalTraits<LO>::one();
139 LO numAggs = aggs->GetNumAggregates();
140 aggSizes = aggs->ComputeAggregateSizes();
142 aggsToIndices = ArrayRCP<LO>(numAggs+LO_ONE,LO_ZERO);
144 for (LO i=0;i<numAggs;++i) {
145 aggsToIndices[i+LO_ONE] = aggsToIndices[i] + aggSizes[i];
148 const RCP<LOMultiVector> vertex2AggId = aggs->GetVertex2AggId();
149 const ArrayRCP<const LO> vertex2AggIdData = vertex2AggId->getData(0);
151 LO numNodes = vertex2AggId->getLocalLength();
152 aggSortedVertices = ArrayRCP<LO>(numNodes,-LO_ONE);
153 std::vector<LO> vertexInsertionIndexByAgg(numNodes,LO_ZERO);
155 for (LO i=0;i<numNodes;++i) {
157 LO aggId = vertex2AggIdData[i];
158 if (aggId<0 || aggId>=numAggs)
continue;
160 aggSortedVertices[aggsToIndices[aggId]+vertexInsertionIndexByAgg[aggId]] = i;
161 vertexInsertionIndexByAgg[aggId]++;
168 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 const SC SCALAR_ONE = Teuchos::ScalarTraits<SC>::one();
172 const SC SCALAR_TWO = SCALAR_ONE + SCALAR_ONE;
174 const LO LO_ZERO = Teuchos::OrdinalTraits<LO>::zero();
175 const LO LO_ONE = Teuchos::OrdinalTraits<LO>::one();
178 const MT MT_ZERO = Teuchos::ScalarTraits<MT>::zero();
179 const MT MT_ONE = Teuchos::ScalarTraits<MT>::one();
180 ParameterList pL = GetParameterList();
182 RCP<const Matrix> AT = A;
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, ALG_REVERSE};
189 if(algostr ==
"forward") {algo = ALG_FORWARD; GetOStream(
Statistics1) <<
"AggregateQuality: Using 'forward' algorithm" << std::endl;}
190 else if(algostr ==
"reverse") {algo = ALG_REVERSE; GetOStream(
Statistics1) <<
"AggregateQuality: Using 'reverse' algorithm" << std::endl;}
195 bool check_symmetry = pL.get<
bool>(
"aggregate qualities: check symmetry");
196 if (check_symmetry) {
198 RCP<MultiVector> x = MultiVectorFactory::Build(A->getMap(), 1,
false);
199 x->Xpetra_randomize();
201 RCP<MultiVector> tmp = MultiVectorFactory::Build(A->getMap(), 1,
false);
203 A->apply(*x, *tmp, Teuchos::NO_TRANS);
204 A->apply(*x, *tmp, Teuchos::TRANS, -SCALAR_ONE, SCALAR_ONE);
206 Array<magnitudeType> tmp_norm(1);
207 tmp->norm2(tmp_norm());
209 Array<magnitudeType> x_norm(1);
210 tmp->norm2(x_norm());
212 if (tmp_norm[0] > 1e-10*x_norm[0]) {
213 std::string transpose_string =
"transpose";
214 RCP<ParameterList> whatever;
217 assert(A->getMap()->isSameAs( *(AT->getMap()) ));
230 ArrayRCP<LO> aggSortedVertices, aggsToIndices, aggSizes;
231 ConvertAggregatesData(aggs, aggSortedVertices, aggsToIndices, aggSizes);
233 LO numAggs = aggs->GetNumAggregates();
237 typedef Teuchos::SerialDenseMatrix<LO,MT> DenseMatrix;
238 typedef Teuchos::SerialDenseVector<LO,MT> DenseVector;
240 ArrayView<const LO> rowIndices;
241 ArrayView<const SC> rowValues;
242 ArrayView<const SC> colValues;
243 Teuchos::LAPACK<LO,MT> myLapack;
246 for (LO aggId=LO_ZERO; aggId<numAggs; ++aggId) {
248 LO aggSize = aggSizes[aggId];
249 DenseMatrix A_aggPart(aggSize, aggSize,
true);
250 DenseVector offDiagonalAbsoluteSums(aggSize,
true);
253 for (LO idx=LO_ZERO; idx<aggSize; ++idx) {
255 LO nodeId = aggSortedVertices[idx+aggsToIndices[aggId]];
256 A->getLocalRowView(nodeId, rowIndices, rowValues);
257 AT->getLocalRowView(nodeId, rowIndices, colValues);
260 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) {
272 if (aggSortedVertices[idx2+aggsToIndices[aggId]] == nodeId2) {
282 if (idxInAgg == -LO_ONE) {
284 offDiagonalAbsoluteSums[idx] += Teuchos::ScalarTraits<SC>::magnitude(val);
288 A_aggPart(idx,idxInAgg) = Teuchos::ScalarTraits<SC>::real(val);
298 DenseMatrix A_aggPartDiagonal(aggSize, aggSize,
true);
299 MT diag_sum = MT_ZERO;
300 for (
int i=0;i<aggSize;++i) {
301 A_aggPartDiagonal(i,i) = Teuchos::ScalarTraits<SC>::real(A_aggPart(i,i));
302 diag_sum += Teuchos::ScalarTraits<SC>::real(A_aggPart(i,i));
305 DenseMatrix ones(aggSize, aggSize,
false);
306 ones.putScalar(MT_ONE);
310 DenseMatrix tmp(aggSize, aggSize,
false);
311 DenseMatrix topMatrix(A_aggPartDiagonal);
313 tmp.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, MT_ONE, ones, A_aggPartDiagonal, MT_ZERO);
314 topMatrix.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -MT_ONE/diag_sum, A_aggPartDiagonal, tmp, MT_ONE);
317 DenseMatrix bottomMatrix(A_aggPart);
318 MT matrixNorm = A_aggPart.normInf();
321 const MT boost = (algo == ALG_FORWARD) ? (-1e4*Teuchos::ScalarTraits<MT>::eps()*matrixNorm) : MT_ZERO;
323 for (
int i=0;i<aggSize;++i){
324 bottomMatrix(i,i) -= offDiagonalAbsoluteSums(i) + boost;
329 DenseVector alpha_real(aggSize,
false);
330 DenseVector alpha_imag(aggSize,
false);
331 DenseVector beta(aggSize,
false);
333 DenseVector workArray(14*(aggSize+1),
false);
335 LO (*ptr2func)(MT*, MT*, MT*);
341 const char compute_flag =
'N';
342 if(algo == ALG_FORWARD) {
344 myLapack.GGES(compute_flag,compute_flag,compute_flag,ptr2func,aggSize,
345 topMatrix.values(),aggSize,bottomMatrix.values(),aggSize,&sdim,
346 alpha_real.values(),alpha_imag.values(),beta.values(),vl,aggSize,
347 vr,aggSize,workArray.values(),workArray.length(),bwork,
349 TEUCHOS_ASSERT(info == LO_ZERO);
351 MT maxEigenVal = MT_ZERO;
352 for (
int i=LO_ZERO;i<aggSize;++i) {
355 maxEigenVal = std::max(maxEigenVal, alpha_real[i]/beta[i]);
358 (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE+MT_ONE)*maxEigenVal;
363 myLapack.GGES(compute_flag,compute_flag,compute_flag,ptr2func,aggSize,
364 bottomMatrix.values(),aggSize,topMatrix.values(),aggSize,&sdim,
365 alpha_real.values(),alpha_imag.values(),beta.values(),vl,aggSize,
366 vr,aggSize,workArray.values(),workArray.length(),bwork,
369 TEUCHOS_ASSERT(info == LO_ZERO);
371 MT minEigenVal = MT_ZERO;
373 for (
int i=LO_ZERO;i<aggSize;++i) {
374 MT ev = alpha_real[i] / beta[i];
375 if(ev > zeroThreshold) {
376 if (minEigenVal == MT_ZERO) minEigenVal = ev;
377 else minEigenVal = std::min(minEigenVal,ev);
380 if(minEigenVal == MT_ZERO) (agg_qualities->getDataNonConst(0))[aggId] = Teuchos::ScalarTraits<MT>::rmax();
381 else (agg_qualities->getDataNonConst(0))[aggId] = (MT_ONE+MT_ONE) / minEigenVal;
386 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
389 ParameterList pL = GetParameterList();
391 magnitudeType good_agg_thresh = Teuchos::as<magnitudeType>(pL.get<
double>(
"aggregate qualities: good aggregate threshold"));
394 ArrayRCP<const MT> data = agg_qualities->getData(0);
399 MT mean_bad_agg = 0.0;
400 MT mean_good_agg = 0.0;
403 for (
size_t i=0;i<agg_qualities->getLocalLength();++i) {
405 if (data[i] > good_agg_thresh) {
407 mean_bad_agg += data[i];
410 mean_good_agg += data[i];
412 worst_agg = std::max(worst_agg, data[i]);
416 if (num_bad_aggs > 0) mean_bad_agg /= num_bad_aggs;
417 mean_good_agg /= agg_qualities->getLocalLength() - num_bad_aggs;
419 if (num_bad_aggs == 0) {
420 GetOStream(
Statistics1) <<
"All aggregates passed the quality measure. Worst aggregate had quality " << worst_agg <<
". Mean aggregate quality " << mean_good_agg <<
"." << std::endl;
422 GetOStream(
Statistics1) << num_bad_aggs <<
" out of " << agg_qualities->getLocalLength() <<
" did not pass the quality measure. Worst aggregate had quality " << worst_agg <<
". " 423 <<
"Mean bad aggregate quality " << mean_bad_agg <<
". Mean good aggregate quality " << mean_good_agg <<
"." << std::endl;
426 if (pL.get<
bool>(
"aggregate qualities: file output")) {
427 std::string filename = pL.get<std::string>(
"aggregate qualities: file base")+
"."+std::to_string(level.
GetLevelID());
428 Xpetra::IO<magnitudeType,LO,GO,Node>::Write(filename, *agg_qualities);
432 const auto n = size_t(agg_qualities->getLocalLength());
437 for (
size_t i=0; i<n; ++i) {
438 tmp.push_back(data[i]);
441 std::sort(tmp.begin(), tmp.end());
443 Teuchos::ArrayView<const double> percents = pL.get<Teuchos::Array<double> >(
"aggregate qualities: percentiles")();
445 GetOStream(
Statistics1) <<
"AGG QUALITY HEADER : | LEVEL | TOTAL |";
446 for (
auto percent : percents) {
447 GetOStream(
Statistics1) << std::fixed << std::setprecision(4) <<100.0*percent <<
"% |";
452 for (
auto percent : percents) {
453 size_t i = size_t(n*percent);
454 i = i < n ? i : n-1u;
456 GetOStream(
Statistics1) << std::fixed <<std::setprecision(4) << tmp[i] <<
" |";
465 #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 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.
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.