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;